158 * <a name=
"elastoplastic.cc-Includefiles"></a>
159 * <h3>Include files</h3>
160 * The
set of include files is not much of a surprise any more at
this time:
163 *
#include <deal.II/base/conditional_ostream.h>
164 *
#include <deal.II/base/parameter_handler.h>
165 *
#include <deal.II/base/utilities.h>
166 *
#include <deal.II/base/index_set.h>
167 *
#include <deal.II/base/quadrature_lib.h>
168 *
#include <deal.II/base/function.h>
169 *
#include <deal.II/base/logstream.h>
170 *
#include <deal.II/base/timer.h>
171 *
#include <deal.II/base/table_handler.h>
173 *
#include <deal.II/lac/vector.h>
174 *
#include <deal.II/lac/full_matrix.h>
175 *
#include <deal.II/lac/sparsity_tools.h>
176 *
#include <deal.II/lac/sparse_matrix.h>
177 *
#include <deal.II/lac/dynamic_sparsity_pattern.h>
178 *
#include <deal.II/lac/block_sparsity_pattern.h>
179 *
#include <deal.II/lac/solver_bicgstab.h>
180 *
#include <deal.II/lac/precondition.h>
181 *
#include <deal.II/lac/affine_constraints.h>
182 *
#include <deal.II/lac/trilinos_sparse_matrix.h>
183 *
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
184 *
#include <deal.II/lac/trilinos_vector.h>
185 *
#include <deal.II/lac/trilinos_precondition.h>
186 *
#include <deal.II/lac/trilinos_solver.h>
187 *
#include <deal.II/lac/sparse_direct.h>
189 *
#include <deal.II/grid/tria.h>
190 *
#include <deal.II/grid/grid_generator.h>
191 *
#include <deal.II/grid/grid_refinement.h>
192 *
#include <deal.II/grid/grid_tools.h>
193 *
#include <deal.II/grid/tria_accessor.h>
194 *
#include <deal.II/grid/tria_iterator.h>
195 *
#include <deal.II/grid/grid_out.h>
196 *
#include <deal.II/grid/manifold_lib.h>
198 *
#include <deal.II/distributed/tria.h>
199 *
#include <deal.II/distributed/grid_refinement.h>
201 *
#include <deal.II/dofs/dof_handler.h>
202 *
#include <deal.II/dofs/dof_accessor.h>
203 *
#include <deal.II/dofs/dof_renumbering.h>
204 *
#include <deal.II/dofs/dof_tools.h>
206 *
#include <deal.II/fe/fe_q.h>
207 *
#include <deal.II/fe/fe_system.h>
208 *
#include <deal.II/fe/fe_values.h>
209 *
#include <deal.II/fe/fe_dgq.h>
210 *
#include <deal.II/fe/fe_tools.h>
212 *
#include <deal.II/numerics/vector_tools.h>
213 *
#include <deal.II/numerics/matrix_tools.h>
214 *
#include <deal.II/numerics/data_out.h>
215 *
#include <deal.II/numerics/error_estimator.h>
216 *
#include <deal.II/numerics/fe_field_function.h>
217 *
#include <deal.II/numerics/solution_transfer.h>
221 * And here the only two
new things among the header files: an include file in
222 * which
symmetric tensors of rank 2 and 4 are implemented, as introduced in
226 *
#include <deal.II/base/symmetric_tensor.h>
230 * And a header that implements filters
for iterators looping over all
231 * cells. We will use
this when selecting only those cells
for output that are
232 * owned by the present process in a %
parallel program:
235 *
#include <deal.II/grid/filtered_iterator.h>
238 *
#include <iostream>
242 * This
final include file provides the <code>mkdir</code> function
243 * that we will use to create a directory
for output files,
if necessary:
246 *
#include <sys/stat.h>
248 *
namespace ElastoPlastic
254 *
const unsigned int n_slices,
260 *
Assert (input.n_levels() == 1,
261 * ExcMessage (
"The input triangulations must be coarse meshes."));
264 *
Assert(result.n_cells()==0, ExcMessage(
"resultin Triangulation need to be empty upon calling extrude_triangulation."));
265 *
Assert(
height>0, ExcMessage(
"The height in extrude_triangulation needs to be positive."));
266 *
Assert(n_slices>=2, ExcMessage(
"The number of slices in extrude_triangulation needs to be at least 2."));
268 *
std::vector<Point<3> > points(n_slices*input.n_used_vertices());
269 *
std::vector<CellData<3> > cells;
270 *
cells.reserve((n_slices-1)*input.n_active_cells());
272 *
for (
unsigned int slice=0; slice<n_slices; ++slice)
274 *
for (
unsigned int i=0; i<input.n_vertices(); ++i)
277 *
if (input.get_used_vertices()[i])
279 *
const Point<2> &v = input.get_vertices()[i];
280 *
points[i+slice*input.n_vertices()](0) = v(0);
281 *
points[i+slice*input.n_vertices()](1) = v(1);
282 *
points[i+slice*input.n_vertices()](2) =
height * slice / (n_slices-1);
288 *
cell = input.begin_active(); cell != input.end(); ++cell)
290 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
293 *
for (
unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
296 *
= cell->vertex_index(v)+slice*input.n_used_vertices();
298 *
= cell->vertex_index(v)+(slice+1)*input.n_used_vertices();
301 *
this_cell.material_id = cell->material_id();
308 *
s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
310 *
cell = input.begin_active(); cell != input.end(); ++cell)
313 *
for (
unsigned int f=0; f<4; ++f)
314 *
if (cell->at_boundary(f))
317 *
bid =
std::max(bid, quad.boundary_id);
318 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
320 *
quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_used_vertices();
321 *
quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_used_vertices();
322 *
quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_used_vertices();
323 *
quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_used_vertices();
324 *
s.boundary_quads.push_back(quad);
330 *
cell = input.begin_active(); cell != input.end(); ++cell)
334 *
quad.vertices[0] = cell->vertex_index(0);
335 *
quad.vertices[1] = cell->vertex_index(1);
336 *
quad.vertices[2] = cell->vertex_index(2);
337 *
quad.vertices[3] = cell->vertex_index(3);
338 *
s.boundary_quads.push_back(quad);
340 *
quad.boundary_id = bid + 2;
341 *
for (
int i=0; i<4; ++i)
342 *
quad.vertices[i] += (n_slices-1)*input.n_used_vertices();
343 *
s.boundary_quads.push_back(quad);
346 *
result.create_triangulation (points,
351 *
namespace Evaluation
363 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
364 * + stress[1][1]*stress[1][1]
365 * - stress[0][0]*stress[1][1]
366 * + 3*stress[0][1]*stress[0][1]);
367 * }
else if (dim == 3)
369 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
370 * + stress[1][1]*stress[1][1]
371 * + stress[2][2]*stress[2][2]
372 * - stress[0][0]*stress[1][1]
373 * - stress[1][1]*stress[2][2]
374 * - stress[0][0]*stress[2][2]
375 * + 3*( stress[0][1]*stress[0][1]
376 * +stress[1][2]*stress[1][2]
377 * +stress[0][2]*stress[0][2]) );
382 * -----------------------------------------------
383 *
"Perforated_strip_tension"
385 *
const double von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
386 * + stress[1][1]*stress[1][1]
387 * - stress[0][0]*stress[1][1]
388 * + 3*stress[0][1]*stress[0][1]);
389 * -----------------------------------------------
391 * plane strain / 3
d case
397 * -----------------------------------------------
409 *
return von_Mises_stress;
414 *
class PointValuesEvaluation
417 *
PointValuesEvaluation (
const Point<dim> &evaluation_point);
425 *
<<
"The evaluation point " << arg1
426 *
<<
" was not found among the vertices of the present grid.");
433 *
PointValuesEvaluation<dim>::
434 *
PointValuesEvaluation (
const Point<dim> &evaluation_point)
436 *
evaluation_point (evaluation_point)
443 *
PointValuesEvaluation<dim>::
448 *
const unsigned int dofs_per_vertex = dof_handler.get_fe().dofs_per_vertex;
450 *
ExcDimensionMismatch (
point_values.size(), dofs_per_vertex));
454 *
cell = dof_handler.begin_active(),
455 *
endc = dof_handler.end();
456 *
bool evaluation_point_found =
false;
457 *
for (; (cell!=endc) && !evaluation_point_found; ++cell)
459 *
if (cell->is_locally_owned() && !evaluation_point_found)
460 *
for (
unsigned int vertex=0;
461 *
vertex<GeometryInfo<dim>::vertices_per_cell;
464 *
if (cell->vertex(vertex).distance (evaluation_point)
466 *
cell->diameter() * 1
e-8)
468 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
470 *
point_values[id] = solution(cell->vertex_dof_index(vertex,
id));
473 *
evaluation_point_found =
true;
480 *
ExcEvaluationPointNotFound(evaluation_point));
489 * <a name=
"elastoplastic.cc-ThecodePointHistorycodeclass"></a>
490 * <h3>The <code>PointHistory</code>
class</h3>
494 * As was mentioned in the introduction, we have to store the old stress in
495 * quadrature
point so that we can compute the residual forces at
this point
496 * during the next time step. This alone would not warrant a structure with
497 * only
one member, but in more complicated applications, we would have to
498 * store more information in quadrature points as well, such as the history
499 * variables of plasticity, etc. In essence, we have to store everything
500 * that affects the present state of the material here, which in plasticity
501 * is determined by the deformation history variables.
505 * We will not give
this class any meaningful functionality beyond being
506 * able to store
data, i.e. there are no constructors, destructors, or other
507 * member
functions. In such cases of `dumb
' classes, we usually opt to
508 * declare them as <code>struct</code> rather than <code>class</code>, to
509 * indicate that they are closer to C-style structures than C++-style
514 * struct PointHistory
516 * SymmetricTensor<2,dim> old_stress;
517 * SymmetricTensor<2,dim> old_strain;
525 * <a name="elastoplastic.cc-ThecodeConstitutiveLawcodeclasstemplate"></a>
526 * <h3>The <code>ConstitutiveLaw</code> class template</h3>
530 * This class provides an interface for a constitutive law, i.e., for the
531 * relationship between strain @f$\varepsilon(\mathbf u)@f$ and stress
532 * @f$\sigma@f$. In this example we are using an elastoplastic material behavior
533 * with linear, isotropic hardening. Such materials are characterized by
534 * Young's modulus @f$E@f$, Poisson
's ratio @f$\nu@f$, the initial yield stress
535 * @f$\sigma_0@f$ and the isotropic hardening parameter @f$\gamma@f$. For @f$\gamma =
536 * 0@f$ we obtain perfect elastoplastic behavior.
540 * As explained in the paper that describes this program, the first Newton
541 * steps are solved with a completely elastic material model to avoid having
542 * to deal with both nonlinearities (plasticity and contact) at once. To this
543 * end, this class has a function <code>set_sigma_0()</code> that we use later
544 * on to simply set @f$\sigma_0@f$ to a very large value -- essentially
545 * guaranteeing that the actual stress will not exceed it, and thereby
546 * producing an elastic material. When we are ready to use a plastic model, we
547 * set @f$\sigma_0@f$ back to its proper value, using the same function. As a
548 * result of this approach, we need to leave <code>sigma_0</code> as the only
549 * non-const member variable of this class.
553 * class ConstitutiveLaw
556 * ConstitutiveLaw (const double E,
558 * const double sigma_0,
559 * const double gamma);
562 * set_sigma_0 (double sigma_zero);
565 * get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
566 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
569 * get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
570 * const std::vector<Tensor<2, dim> > &point_hessian,
571 * Tensor<5, dim> &stress_strain_tensor_grad) const;
574 * get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor,
575 * SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
576 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
579 * const double kappa;
582 * const double gamma;
584 * const SymmetricTensor<4, dim> stress_strain_tensor_kappa;
585 * const SymmetricTensor<4, dim> stress_strain_tensor_mu;
590 * The constructor of the ConstitutiveLaw class sets the required material
591 * parameter for our deformable body. Material parameters for elastic
592 * isotropic media can be defined in a variety of ways, such as the pair @f$E,
593 * \nu@f$ (elastic modulus and Poisson's number),
using the Lame parameters
594 * @f$\lambda,mu@f$ or several other commonly used conventions. Here, the
595 * constructor takes a description of material parameters in the form of
596 * @f$E,\nu@f$, but since
this turns out to these are not the coefficients that
597 * appear in the equations of the plastic projector, we immediately convert
598 * them into the more suitable
set @f$\kappa,\mu@f$ of bulk and shear moduli. In
599 * addition, the constructor takes @f$\sigma_0@f$ (the yield stress absent any
600 * plastic strain) and @f$\gamma@f$ (the hardening parameter) as arguments. In
601 *
this constructor, we also compute the two principal components of the
602 * stress-strain relation and its linearization.
606 *
ConstitutiveLaw<dim>::ConstitutiveLaw (
double E,
613 * --------------------
615 * kappa (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (3 * (1 - 2 * (nu / (1+nu))))),
616 *
mu (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))),
617 * --------------------
618 * 3
d and plane strain
621 *
kappa (E / (3 * (1 - 2 * nu))),
622 *
mu (E / (2 * (1 + nu))),
625 * --------------------
630 *
stress_strain_tensor_kappa (kappa
632 *
unit_symmetric_tensor<dim>())),
633 *
stress_strain_tensor_mu (2 * mu
634 *
* (identity_tensor<dim>()
636 *
unit_symmetric_tensor<dim>()) / 3.0))
642 *
ConstitutiveLaw<dim>::set_sigma_0 (
double sigma_zero)
644 *
sigma_0 = sigma_zero;
651 * <a name=
"elastoplastic.cc-ConstitutiveLawget_stress_strain_tensor"></a>
652 * <h4>ConstitutiveLaw::get_stress_strain_tensor</h4>
656 * This is the principal
component of the constitutive law. It projects the
657 * deviatoric part of the stresses in a quadrature
point back to the yield
658 * stress (i.e., the original yield stress @f$\sigma_0@f$ plus the term that
659 * describes linear isotropic hardening). We need
this function to calculate
660 * the nonlinear residual in PlasticityContactProblem::residual_nl_system. The
661 * computations follow the formulas laid out in the introduction.
665 * The function returns whether the quadrature
point is plastic to allow
for
666 * some statistics
downstream on how many of the quadrature points are
667 * plastic and how many are elastic.
672 *
ConstitutiveLaw<dim>::
677 *
stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
683 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
686 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
688 *
stress_strain_tensor = stress_strain_tensor_mu;
689 *
if (von_Mises_stress > sigma_0)
691 *
const double beta = sigma_0 / von_Mises_stress;
692 *
stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
695 *
stress_strain_tensor += stress_strain_tensor_kappa;
697 *
return (von_Mises_stress > sigma_0);
703 *
ConstitutiveLaw<dim>::
709 *
stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
712 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
714 *
if (von_Mises_stress > sigma_0)
717 *
const double deviator_strain_tensor_norm = deviator_strain_tensor.norm();
718 *
const double multiplier = -(1-
gamma)*sigma_0/(2*mu*
std::pow(deviator_strain_tensor_norm,3));
721 *
multiplier_vector = 0;
723 *
for (
unsigned int i=0; i!=dim; ++i)
724 *
for (
unsigned int m=0; m!=dim; ++m)
725 *
for (
unsigned int n=0; n!=dim; ++n)
727 *
multiplier_vector(i) += deviator_strain_tensor[m][n] *
728 *
( 0.5*( point_hessian[m][n][i] + point_hessian[n][m][i] )
729 *
+ ( m==n && dim==2 ? -1/dim*(point_hessian[0][0][i]
730 *
+ point_hessian[1][1][i]) : 0 )
731 *
+ ( m==n && dim==3 ? -1/dim*(point_hessian[0][0][i]
732 *
+ point_hessian[1][1][i]
733 *
+ point_hessian[2][2][i]) : 0 ) );
738 * -----------------------------------------------
739 *
"Perforated_strip_tension"
742 * -----------------------------------------------
744 * plane strain / 3
d case
747 *
const double VM_factor =
std::sqrt(1.5);
750 * -----------------------------------------------
756 *
for (
unsigned int i=0; i!=dim; ++i)
757 *
for (
unsigned int j=0; j!=dim; ++j)
758 *
for (
unsigned int k=0; k!=dim; ++k)
759 *
for (
unsigned int l=0;
l!=dim; ++
l)
760 *
for (
unsigned int m=0; m!=dim; ++m)
762 *
stress_strain_tensor_grad[i][j][k][
l][m] = 1/VM_factor
764 *
* stress_strain_tensor_mu[i][j][k][
l]
765 *
* multiplier_vector(m);
771 *
stress_strain_tensor_grad = 0;
774 *
return (von_Mises_stress > sigma_0);
781 * <a name=
"elastoplastic.cc-ConstitutiveLawget_linearized_stress_strain_tensors"></a>
782 * <h4>ConstitutiveLaw::get_linearized_stress_strain_tensors</h4>
786 * This function returns the linearized stress strain tensor, linearized
787 * around the solution @f$u^{i-1}@f$ of the previous Newton step @f$i-1@f$. The
788 * parameter <code>strain_tensor</code> (commonly denoted
789 * @f$\varepsilon(u^{i-1})@f$) must be passed as an argument, and serves as the
790 * linearization
point. The function returns the derivative of the nonlinear
791 * constitutive law in the variable stress_strain_tensor, as well as the
792 * stress-strain tensor of the linearized problem in
793 * stress_strain_tensor_linearized. See
794 * PlasticityContactProblem::assemble_nl_system where
this function is used.
799 *
ConstitutiveLaw<dim>::
805 *
stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
808 *
stress_strain_tensor = stress_strain_tensor_mu;
809 *
stress_strain_tensor_linearized = stress_strain_tensor_mu;
812 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
813 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
815 *
if (von_Mises_stress > sigma_0)
817 *
const double beta = sigma_0 / von_Mises_stress;
818 *
stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
819 *
stress_strain_tensor_linearized *= (
gamma + (1 -
gamma) * beta);
820 *
deviator_stress_tensor /= deviator_stress_tensor_norm;
821 *
stress_strain_tensor_linearized -= (1 -
gamma) * beta * 2 * mu
823 *
deviator_stress_tensor);
826 *
stress_strain_tensor += stress_strain_tensor_kappa;
827 *
stress_strain_tensor_linearized += stress_strain_tensor_kappa;
832 * Finally, below we will need a function that computes the rotation
matrix
833 * induced by a displacement at a given
point. In fact, of course, the
834 * displacement at a single
point only has a direction and a magnitude, it
835 * is the change in direction and magnitude that induces rotations. In
837 * displacement, or, more specifically, from the curl.
841 * The formulas by which the rotation matrices are determined are a little
842 * awkward, especially in 3
d. For 2
d, there is a simpler way, so we
843 * implement
this function twice, once
for 2
d and once
for 3
d, so that we
844 * can compile and use the program in both space dimensions
if so desired --
845 * after all, deal.II is all about dimension
independent programming and
846 * reuse of algorithm thoroughly tested with cheap computations in 2
d,
for
847 * the more expensive computations in 3
d. Here is
one case, where we have to
848 * implement different algorithms
for 2
d and 3
d, but then can write the rest
849 * of the program in a way that is
independent of the space dimension.
853 * So, without further ado to the 2
d implementation:
857 *
get_rotation_matrix (
const std::vector<
Tensor<1,2> > &grad_u)
861 * First, compute the curl of the velocity field from the
gradients. Note
862 * that we are in 2
d, so the rotation is a
scalar:
865 *
const double curl = (grad_u[1][0] - grad_u[0][1]);
869 * From
this, compute the
angle of rotation:
876 * And from
this, build the antisymmetric rotation
matrix:
879 *
const double t[2][2] = {{
cos(angle),
sin(angle) },
880 *
{-
sin(angle),
cos(angle) }
888 * The 3
d case is a little more contrived:
892 *
get_rotation_matrix (
const std::vector<
Tensor<1,3> > &grad_u)
896 * Again
first compute the curl of the velocity field. This time, it is a
900 *
const Point<3> curl (grad_u[2][1] - grad_u[1][2],
901 *
grad_u[0][2] - grad_u[2][0],
902 *
grad_u[1][0] - grad_u[0][1]);
906 * From
this vector,
using its magnitude, compute the tangent of the
angle
907 * of rotation, and from it the actual
angle:
910 *
const double tan_angle =
std::sqrt(curl*curl);
915 * Now, here
's one problem: if the angle of rotation is too small, that
916 * means that there is no rotation going on (for example a translational
917 * motion). In that case, the rotation matrix is the identity matrix.
921 * The reason why we stress that is that in this case we have that
922 * <code>tan_angle==0</code>. Further down, we need to divide by that
923 * number in the computation of the axis of rotation, and we would get
924 * into trouble when dividing doing so. Therefore, let's shortcut
this and
925 * simply
return the identity
matrix if the
angle of rotation is really
931 *
static const double rotation[3][3]
932 *
= {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } };
939 * Otherwise compute the real rotation
matrix. The algorithm
for this is
940 * not exactly obvious, but can be found in a number of books,
941 * particularly on computer games where rotation is a very frequent
942 * operation. Online, you can find a description at
944 * the signs as here) at
950 *
const double t = 1-c;
952 *
const Point<3> axis = curl/tan_angle;
953 *
const double rotation[3][3]
955 *
t *axis[0] *axis[0]+c,
956 *
t *axis[0] *axis[1]+s *axis[2],
957 *
t *axis[0] *axis[2]-s *axis[1]
960 *
t *axis[0] *axis[1]-s *axis[2],
961 *
t *axis[1] *axis[1]+c,
962 *
t *axis[1] *axis[2]+s *axis[0]
965 *
t *axis[0] *axis[2]+s *axis[1],
966 *
t *axis[1] *axis[1]-s *axis[0],
967 *
t *axis[2] *axis[2]+c
976 * <h3>Equation
data: Body forces, boundary forces,
977 * incremental boundary
values</h3>
981 * The following should be relatively standard. We need classes
for
982 * the boundary forcing term (which we here choose to be zero)
983 * and incremental boundary
values.
986 *
namespace EquationData
1232 * ----------------------------- TimoshenkoBeam ---------------------------------------
1435 * ------------------------- Thick_tube_internal_pressure ----------------------------------
1604 * ------------------------- Perforated_strip_tension ----------------------------------
1774 * ------------------------- Cantiliver_beam_3d ----------------------------------
1777 *
template <
int dim>
1778 *
class IncrementalBoundaryForce :
public Function<dim>
1781 *
IncrementalBoundaryForce (
const double present_time,
1782 *
const double end_time);
1794 *
const double present_time,
1800 *
template <
int dim>
1801 *
IncrementalBoundaryForce<dim>::
1802 *
IncrementalBoundaryForce (
const double present_time,
1803 *
const double end_time)
1806 *
present_time (present_time),
1807 *
end_time (end_time),
1812 *
template <
int dim>
1814 *
IncrementalBoundaryForce<dim>::vector_value (
const Point<dim> &p,
1819 *
ExcDimensionMismatch (
values.size(), dim));
1825 * pressure should be imposed on the top surface, y =
height
1836 * compute the fraction of imposed force
1839 *
const double frac = present_time/end_time;
1844 *
template <
int dim>
1846 *
IncrementalBoundaryForce<dim>::
1847 *
vector_value_list (
const std::vector<
Point<dim> > &points,
1850 *
const unsigned int n_points = points.size();
1852 *
Assert (value_list.size() == n_points,
1853 *
ExcDimensionMismatch (value_list.size(), n_points));
1855 *
for (
unsigned int p=0; p<n_points; ++p)
1856 *
IncrementalBoundaryForce<dim>::vector_value (points[p], value_list[p]);
1860 *
template <
int dim>
1864 *
BodyForce () :
Functions::ZeroFunction<dim> (dim) {}
1868 *
template <
int dim>
1869 *
class IncrementalBoundaryValues :
public Function<dim>
1872 *
IncrementalBoundaryValues (
const double present_time,
1873 *
const double end_time);
1886 *
const double present_time,
1891 *
template <
int dim>
1892 *
IncrementalBoundaryValues<dim>::
1893 *
IncrementalBoundaryValues (
const double present_time,
1894 *
const double end_time)
1897 *
present_time (present_time),
1898 *
end_time (end_time)
1902 *
template <
int dim>
1904 *
IncrementalBoundaryValues<dim>::
1909 *
ExcDimensionMismatch (
values.size(), dim));
1916 *
template <
int dim>
1918 *
IncrementalBoundaryValues<dim>::
1919 *
vector_value_list (
const std::vector<
Point<dim> > &points,
1922 *
const unsigned int n_points = points.size();
1924 *
Assert (value_list.size() == n_points,
1925 *
ExcDimensionMismatch (value_list.size(), n_points));
1927 *
for (
unsigned int p=0; p<n_points; ++p)
1928 *
IncrementalBoundaryValues<dim>::vector_value (points[p], value_list[p]);
1933 * -------------------------------------------------------------------------------
1939 *
namespace DualFunctional
1942 *
template <
int dim>
1950 *
const ConstitutiveLaw<dim> &constitutive_law,
1956 *
template <
int dim>
1957 *
class PointValuesEvaluation :
public DualFunctionalBase<dim>
1960 *
PointValuesEvaluation (
const Point<dim> &evaluation_point);
1966 *
const ConstitutiveLaw<dim> &constitutive_law,
1972 *
<<
"The evaluation point " << arg1
1973 *
<<
" was not found among the vertices of the present grid.");
1980 *
template <
int dim>
1981 *
PointValuesEvaluation<dim>::
1982 *
PointValuesEvaluation (
const Point<dim> &evaluation_point)
1984 *
evaluation_point (evaluation_point)
1988 *
template <
int dim>
1990 *
PointValuesEvaluation<dim>::
1993 *
const ConstitutiveLaw<dim> &,
1997 *
rhs_dual.reinit (dof_handler_dual.n_dofs());
1998 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2001 *
cell_dual = dof_handler_dual.begin_active(),
2002 *
endc_dual = dof_handler_dual.end();
2003 *
for (; cell_dual!=endc_dual; ++cell_dual)
2004 *
for (
unsigned int vertex=0;
2005 *
vertex<GeometryInfo<dim>::vertices_per_cell;
2007 *
if (cell_dual->vertex(vertex).distance(evaluation_point)
2010 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
2012 *
rhs_dual(cell_dual->vertex_dof_index(vertex,
id)) = 1;
2017 *
AssertThrow (
false, ExcEvaluationPointNotFound(evaluation_point));
2021 *
template <
int dim>
2022 *
class PointXDerivativesEvaluation :
public DualFunctionalBase<dim>
2025 *
PointXDerivativesEvaluation (
const Point<dim> &evaluation_point);
2031 *
const ConstitutiveLaw<dim> &constitutive_law,
2037 *
<<
"The evaluation point " << arg1
2038 *
<<
" was not found among the vertices of the present grid.");
2045 *
template <
int dim>
2046 *
PointXDerivativesEvaluation<dim>::
2047 *
PointXDerivativesEvaluation (
const Point<dim> &evaluation_point)
2049 *
evaluation_point (evaluation_point)
2053 *
template <
int dim>
2055 *
PointXDerivativesEvaluation<dim>::
2058 *
const ConstitutiveLaw<dim> &,
2062 *
rhs_dual.reinit (dof_handler_dual.n_dofs());
2065 *
FEValues<dim> fe_values (dof_handler_dual.get_fe(), quadrature,
2069 *
const unsigned int n_q_points = fe_values.n_quadrature_points;
2070 *
Assert ( n_q_points==quadrature.size() , ExcInternalError() );
2071 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2074 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2076 *
double total_volume = 0;
2079 *
cell = dof_handler_dual.begin_active(),
2080 *
endc = dof_handler_dual.end();
2081 *
for (; cell!=endc; ++cell)
2082 *
if (cell->center().distance(evaluation_point) <=
2085 *
fe_values.reinit (cell);
2088 *
for (
unsigned int q=0; q<n_q_points; ++q)
2090 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2092 *
cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
2093 *
fe_values.JxW (q);
2096 *
total_volume += fe_values.JxW (q);
2099 *
cell->get_dof_indices (local_dof_indices);
2100 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2102 *
rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2107 *
ExcEvaluationPointNotFound(evaluation_point));
2109 *
rhs_dual *= 1./total_volume;
2114 *
template <
int dim>
2115 *
class MeanDisplacementFace :
public DualFunctionalBase<dim>
2118 *
MeanDisplacementFace (
const unsigned int face_id,
2119 *
const std::vector<bool> comp_mask);
2125 *
const ConstitutiveLaw<dim> &constitutive_law,
2130 *
const unsigned int face_id;
2131 *
const std::vector<bool> comp_mask;
2135 *
template <
int dim>
2136 *
MeanDisplacementFace<dim>::
2137 *
MeanDisplacementFace (
const unsigned int face_id,
2138 *
const std::vector<bool> comp_mask )
2140 *
face_id (face_id),
2141 *
comp_mask (comp_mask)
2144 *
ExcDimensionMismatch (comp_mask.size(), dim) );
2148 *
template <
int dim>
2150 *
MeanDisplacementFace<dim>::
2153 *
const ConstitutiveLaw<dim> &,
2159 *
rhs_dual.reinit (dof_handler_dual.n_dofs());
2161 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2165 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2166 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2167 *
const unsigned int n_face_q_points = face_quadrature.size();
2170 *
ExcDimensionMismatch (dofs_per_vertex, dim) );
2172 *
std::vector<unsigned int> comp_vector(dofs_per_vertex);
2173 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2177 *
comp_vector[i] = 1;
2183 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2187 * bound_size :
size of the boundary, in 2
d is the length
2188 * and in the 3
d case, area
2191 *
double bound_size = 0.;
2194 *
cell = dof_handler_dual.begin_active(),
2195 *
endc = dof_handler_dual.end();
2196 *
bool evaluation_face_found =
false;
2197 *
for (; cell!=endc; ++cell)
2200 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2202 *
if (cell->face(face)->at_boundary()
2204 *
cell->face(face)->boundary_id() == face_id)
2206 *
if (!evaluation_face_found)
2208 *
evaluation_face_found =
true;
2210 *
fe_face_values.reinit (cell, face);
2212 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2214 *
bound_size += fe_face_values.JxW(q_point);
2216 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2218 *
const unsigned int
2219 *
component_i = dof_handler_dual.get_fe().system_to_component_index(i).first;
2221 *
cell_rhs(i) += (fe_face_values.shape_value(i,q_point) *
2222 *
comp_vector[component_i] *
2223 *
fe_face_values.JxW(q_point));
2231 *
cell->get_dof_indices (local_dof_indices);
2232 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2234 *
rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2239 *
AssertThrow(evaluation_face_found, ExcInternalError());
2241 *
rhs_dual /= bound_size;
2246 *
template <
int dim>
2247 *
class MeanStressFace :
public DualFunctionalBase<dim>
2250 *
MeanStressFace (
const unsigned int face_id,
2251 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2257 *
const ConstitutiveLaw<dim> &constitutive_law,
2262 *
const unsigned int face_id;
2263 *
const std::vector<std::vector<unsigned int> > comp_stress;
2267 *
template <
int dim>
2268 *
MeanStressFace<dim>::
2269 *
MeanStressFace (
const unsigned int face_id,
2270 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2272 *
face_id (face_id),
2273 *
comp_stress (comp_stress)
2276 *
ExcDimensionMismatch (comp_stress.size(), dim) );
2280 *
template <
int dim>
2282 *
MeanStressFace<dim>::
2285 *
const ConstitutiveLaw<dim> &constitutive_law,
2291 *
rhs_dual.reinit (dof_handler_dual.n_dofs());
2293 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2297 *
FEFaceValues<dim> fe_face_values_dual (dof_handler_dual.get_fe(), face_quadrature,
2300 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2301 *
const unsigned int n_face_q_points = face_quadrature.size();
2303 *
std::vector<SymmetricTensor<2, dim> > strain_tensor(n_face_q_points);
2308 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2312 * bound_size :
size of the boundary, in 2
d is the length
2313 * and in the 3
d case, area
2316 *
double bound_size = 0.;
2318 *
bool evaluation_face_found =
false;
2321 *
cell_dual = dof_handler_dual.begin_active(),
2322 *
endc_dual = dof_handler_dual.end(),
2323 *
cell = dof_handler.begin_active();
2327 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2330 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2332 *
if (cell_dual->face(face)->at_boundary()
2334 *
cell_dual->face(face)->boundary_id() == face_id)
2336 *
if (!evaluation_face_found)
2338 *
evaluation_face_found =
true;
2341 *
fe_face_values.reinit (cell, face);
2342 *
fe_face_values_dual.reinit (cell_dual, face);
2344 *
fe_face_values[displacement].get_function_symmetric_gradients(solution,
2347 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2349 *
bound_size += fe_face_values_dual.JxW(q_point);
2351 *
constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2352 *
stress_strain_tensor);
2354 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2357 *
stress_phi_i = stress_strain_tensor
2358 *
* fe_face_values_dual[displacement].symmetric_gradient(i, q_point);
2360 *
for (
unsigned int k=0; k!=dim; ++k)
2362 *
for (
unsigned int l=0;
l!=dim; ++
l)
2364 *
if ( comp_stress[k][l] == 1 )
2366 *
cell_rhs(i) += stress_phi_i[k][
l]
2368 *
fe_face_values_dual.JxW(q_point);
2381 *
cell_dual->get_dof_indices (local_dof_indices);
2382 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2384 *
rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2389 *
AssertThrow(evaluation_face_found, ExcInternalError());
2391 *
rhs_dual /= bound_size;
2396 *
template <
int dim>
2397 *
class MeanStressDomain :
public DualFunctionalBase<dim>
2400 *
MeanStressDomain (
const std::string &base_mesh,
2401 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2407 *
const ConstitutiveLaw<dim> &constitutive_law,
2412 *
const std::string base_mesh;
2413 *
const std::vector<std::vector<unsigned int> > comp_stress;
2417 *
template <
int dim>
2418 *
MeanStressDomain<dim>::
2419 *
MeanStressDomain (
const std::string &base_mesh,
2420 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2422 *
base_mesh (base_mesh),
2423 *
comp_stress (comp_stress)
2426 *
ExcDimensionMismatch (comp_stress.size(), dim) );
2430 *
template <
int dim>
2432 *
MeanStressDomain<dim>::
2435 *
const ConstitutiveLaw<dim> &constitutive_law,
2439 *
AssertThrow (base_mesh ==
"Cantiliver_beam_3d", ExcNotImplemented());
2444 * Mean stress at the specified domain is of interest.
2445 * The interest domains are located on the bottom and top of the flanges
2446 * close to the clamped face, z = 0
2447 * top domain:
height/2 - thickness_flange <= y <=
height/2
2448 * 0 <= z <= 2 * thickness_flange
2449 * bottom domain: -
height/2 <= y <= -
height/2 + thickness_flange
2450 * 0 <= z <= 2 * thickness_flange
2456 *
const double height = 200
e-3,
2457 *
thickness_flange = 10
e-3;
2459 *
rhs_dual.reinit (dof_handler_dual.n_dofs());
2461 *
const QGauss<dim> quadrature_formula(dof_handler_dual.get_fe().tensor_degree()+1);
2463 *
FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
2465 *
FEValues<dim> fe_values_dual (dof_handler_dual.get_fe(), quadrature_formula,
2468 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2469 *
const unsigned int n_q_points = quadrature_formula.size();
2471 *
std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
2476 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2480 * domain_size :
size of the interested domain, in 2
d is the area
2484 *
double domain_size = 0.;
2486 *
bool evaluation_domain_found =
false;
2489 *
cell_dual = dof_handler_dual.begin_active(),
2490 *
endc_dual = dof_handler_dual.end(),
2491 *
cell = dof_handler.begin_active();
2495 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2497 *
const double y = cell->center()[1],
2498 *
z = cell->center()[2];
2501 * top domain:
height/2 - thickness_flange <= y <=
height/2
2502 * 0 <= z <= 2 * thickness_flange
2503 * bottom domain: -
height/2 <= y <= -
height/2 + thickness_flange
2504 * 0 <= z <= 2 * thickness_flange
2507 *
if ( ((z > 0) && (z < 2*thickness_flange)) &&
2508 *
( ((y >
height/2 - thickness_flange) && (y <
height/2)) ||
2509 *
((y > -
height/2) && (y < -
height/2 + thickness_flange)) ) )
2513 *
if (!evaluation_domain_found)
2515 *
evaluation_domain_found =
true;
2518 *
fe_values.reinit(cell);
2519 *
fe_values_dual.reinit(cell_dual);
2521 *
fe_values[displacement].get_function_symmetric_gradients(solution,
2524 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
2526 *
domain_size += fe_values_dual.JxW(q_point);
2528 *
constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2529 *
stress_strain_tensor);
2531 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2534 *
stress_phi_i = stress_strain_tensor
2535 *
* fe_values_dual[displacement].symmetric_gradient(i, q_point);
2537 *
for (
unsigned int k=0; k!=dim; ++k)
2539 *
for (
unsigned int l=0;
l!=dim; ++
l)
2541 *
if ( comp_stress[k][l] == 1 )
2543 *
cell_rhs(i) += stress_phi_i[k][
l]
2545 *
fe_values_dual.JxW(q_point);
2557 *
cell_dual->get_dof_indices (local_dof_indices);
2558 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2560 *
rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2565 *
AssertThrow(evaluation_domain_found, ExcInternalError());
2567 *
rhs_dual /= domain_size;
2572 *
template <
int dim>
2573 *
class MeanStrainEnergyFace :
public DualFunctionalBase<dim>
2576 *
MeanStrainEnergyFace (
const unsigned int face_id,
2580 *
void assemble_rhs_nonlinear (
const DoFHandler<dim> &primal_dof_handler,
2586 *
const unsigned int face_id;
2592 *
template <
int dim>
2593 *
MeanStrainEnergyFace<dim>::
2594 *
MeanStrainEnergyFace (
const unsigned int face_id,
2598 *
face_id (face_id),
2599 *
lambda_function (&lambda_function),
2600 *
mu_function (&mu_function)
2604 *
template <
int dim>
2606 *
MeanStrainEnergyFace<dim>::
2614 * Assemble right hand side of the dual problem when the quantity of interest is
2615 * a nonlinear functional. In
this case, the QoI should be linearized which depends
2616 * on the solution of the primal problem.
2617 * The extractor of the linearized QoI functional is the
gradient of the the original
2618 * QoI functional with the primal solution
values.
2626 *
rhs.reinit (dof_handler.n_dofs());
2628 *
const QGauss<dim-1> face_quadrature(dof_handler.get_fe().tensor_degree()+1);
2629 *
FEFaceValues<dim> primal_fe_face_values (primal_dof_handler.get_fe(), face_quadrature,
2637 *
const unsigned int dofs_per_vertex = primal_dof_handler.get_fe().dofs_per_vertex;
2638 *
const unsigned int n_face_q_points = face_quadrature.size();
2639 *
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
2642 *
ExcDimensionMismatch (dofs_per_vertex, dim) );
2644 *
std::vector< std::vector< Tensor<1,dim> > > primal_solution_gradients;
2645 *
primal_solution_gradients.resize(n_face_q_points);
2647 *
std::vector<std::vector<Tensor<2,dim> > > primal_solution_hessians;
2648 *
primal_solution_hessians.resize (n_face_q_points);
2650 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
2652 *
primal_solution_gradients[i].resize (dofs_per_vertex);
2653 *
primal_solution_hessians[i].resize (dofs_per_vertex);
2656 *
std::vector<double> lambda_values (n_face_q_points);
2657 *
std::vector<double> mu_values (n_face_q_points);
2661 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2665 * bound_size :
size of the boundary, in 2
d is the length
2666 * and in the 3
d case, area
2669 *
double bound_size = 0.;
2671 *
bool evaluation_face_found =
false;
2674 *
primal_cell = primal_dof_handler.begin_active(),
2675 *
primal_endc = primal_dof_handler.end();
2678 *
cell = dof_handler.begin_active(),
2679 *
endc = dof_handler.end();
2681 *
for (; cell!=endc; ++cell, ++primal_cell)
2684 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2686 *
if (cell->face(face)->at_boundary()
2688 *
cell->face(face)->boundary_id() == face_id)
2690 *
if (!evaluation_face_found)
2692 *
evaluation_face_found =
true;
2694 *
primal_fe_face_values.reinit (primal_cell, face);
2696 *
primal_fe_face_values.get_function_gradients (primal_solution,
2697 *
primal_solution_gradients);
2699 *
primal_fe_face_values.get_function_hessians (primal_solution,
2700 *
primal_solution_hessians);
2702 *
lambda_function->value_list (primal_fe_face_values.get_quadrature_points(), lambda_values);
2703 *
mu_function->value_list (primal_fe_face_values.get_quadrature_points(), mu_values);
2705 *
fe_face_values.reinit (cell, face);
2707 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2709 *
bound_size += primal_fe_face_values.JxW(q_point);
2711 *
for (
unsigned int m=0; m<dofs_per_cell; ++m)
2713 *
const unsigned int
2714 *
component_m = dof_handler.get_fe().system_to_component_index(m).first;
2716 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2718 *
for (
unsigned int j=0; j!=dofs_per_vertex; ++j)
2720 *
cell_rhs(m) += fe_face_values.shape_value(m,q_point) *
2722 *
lambda_values[q_point] *
2724 *
primal_solution_hessians[q_point][i][i][component_m] * primal_solution_gradients[q_point][j][j]
2726 *
primal_solution_gradients[q_point][i][i] * primal_solution_hessians[q_point][j][j][component_m]
2729 *
mu_values[q_point] *
2731 *
2*primal_solution_hessians[q_point][j][i][component_m] * primal_solution_gradients[q_point][j][i]
2733 *
primal_solution_hessians[q_point][i][j][component_m] * primal_solution_gradients[q_point][j][i]
2735 *
primal_solution_gradients[q_point][i][j] * primal_solution_hessians[q_point][j][i][component_m]
2738 *
primal_fe_face_values.JxW(q_point);
2751 *
cell->get_dof_indices (local_dof_indices);
2752 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2754 *
rhs(local_dof_indices[i]) += cell_rhs(i);
2759 *
AssertThrow(evaluation_face_found, ExcInternalError());
2761 *
rhs *= 1./(2*bound_size);
2774 *
template <
int dim>
2781 *
const ConstitutiveLaw<dim> &constitutive_law,
2782 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2783 *
const unsigned int ×tep_no,
2784 *
const std::string &output_dir,
2785 *
const std::string &base_mesh,
2786 *
const double &present_time,
2787 *
const double &end_time);
2789 *
void compute_error_DWR (
Vector<float> &estimated_error_per_cell);
2794 *
void setup_system ();
2795 *
void compute_dirichlet_constraints ();
2796 *
void assemble_matrix ();
2797 *
void assemble_rhs ();
2799 *
void output_results ();
2805 *
const unsigned int fe_degree;
2808 *
const unsigned int fe_degree_dual;
2813 *
const QGauss<dim - 1> face_quadrature_formula;
2823 *
const ConstitutiveLaw<dim> constitutive_law;
2828 *
unsigned int timestep_no;
2829 *
std::string output_dir;
2830 *
const std::string base_mesh;
2831 *
double present_time;
2841 *
const ConstitutiveLaw<dim> &constitutive_law,
2842 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2843 *
const unsigned int ×tep_no,
2844 *
const std::string &output_dir,
2845 *
const std::string &base_mesh,
2846 *
const double &present_time,
2847 *
const double &end_time)
2850 *
dof_handler (triangulation),
2851 *
solution(solution),
2852 *
fe_degree(fe.tensor_degree()),
2853 *
fe_degree_dual(fe_degree + 1),
2854 *
fe_dual(
FE_Q<dim>(fe_degree_dual), dim),
2855 *
dof_handler_dual (triangulation),
2856 *
quadrature_formula (fe_degree_dual + 1),
2857 *
face_quadrature_formula (fe_degree_dual + 1),
2858 *
constitutive_law (constitutive_law),
2859 *
triangulation (&triangulation),
2860 *
dual_functional (&dual_functional),
2861 *
timestep_no (timestep_no),
2862 *
output_dir (output_dir),
2863 *
base_mesh (base_mesh),
2864 *
present_time (present_time),
2865 *
end_time (end_time)
2870 *
DualSolver<dim>::~DualSolver()
2872 *
dof_handler_dual.
clear ();
2877 *
void DualSolver<dim>::setup_system()
2879 *
dof_handler.distribute_dofs(fe);
2881 *
dof_handler_dual.distribute_dofs (fe_dual);
2882 *
std::cout <<
" Number of degrees of freedom in dual problem: "
2883 *
<< dof_handler_dual.n_dofs()
2886 *
constraints_hanging_nodes_dual.clear ();
2888 *
constraints_hanging_nodes_dual);
2889 *
constraints_hanging_nodes_dual.close ();
2891 *
compute_dirichlet_constraints();
2893 *
sparsity_pattern_dual.reinit (dof_handler_dual.n_dofs(),
2894 *
dof_handler_dual.n_dofs(),
2895 *
dof_handler_dual.max_couplings_between_dofs());
2900 * constraints_hanging_nodes_dual.condense (sparsity_pattern_dual);
2903 *
constraints_dirichlet_and_hanging_nodes_dual.condense (sparsity_pattern_dual);
2905 *
sparsity_pattern_dual.compress();
2907 *
system_matrix_dual.reinit (sparsity_pattern_dual);
2909 *
solution_dual.reinit (dof_handler_dual.n_dofs());
2910 *
system_rhs_dual.reinit (dof_handler_dual.n_dofs());
2915 *
void DualSolver<dim>::compute_dirichlet_constraints()
2917 *
constraints_dirichlet_and_hanging_nodes_dual.clear ();
2918 *
constraints_dirichlet_and_hanging_nodes_dual.merge(constraints_hanging_nodes_dual);
2920 *
std::vector<bool> component_mask(dim);
2922 *
if (base_mesh ==
"Timoshenko beam")
2926 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2927 *
constraints_dirichlet_and_hanging_nodes_dual,
2930 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
2934 * the boundary x = 0
2937 *
component_mask[0] =
true;
2938 *
component_mask[1] =
false;
2941 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2942 *
constraints_dirichlet_and_hanging_nodes_dual,
2946 * the boundary y = 0
2949 *
component_mask[0] =
false;
2950 *
component_mask[1] =
true;
2953 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2954 *
constraints_dirichlet_and_hanging_nodes_dual,
2957 *
else if (base_mesh ==
"Perforated_strip_tension")
2961 * the boundary x = 0
2964 *
component_mask[0] =
true;
2965 *
component_mask[1] =
false;
2966 *
component_mask[2] =
false;
2969 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2970 *
constraints_dirichlet_and_hanging_nodes_dual,
2974 * the boundary y = 0
2977 *
component_mask[0] =
false;
2978 *
component_mask[1] =
true;
2979 *
component_mask[2] =
false;
2982 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2983 *
constraints_dirichlet_and_hanging_nodes_dual,
2987 * the boundary y = imposed incremental displacement
2990 *
component_mask[0] =
false;
2991 *
component_mask[1] =
true;
2992 *
component_mask[2] =
false;
2995 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2996 *
constraints_dirichlet_and_hanging_nodes_dual,
2999 *
else if (base_mesh ==
"Cantiliver_beam_3d")
3003 * the boundary x = y = z = 0
3006 *
component_mask[0] =
true;
3007 *
component_mask[1] =
true;
3008 *
component_mask[2] =
true;
3011 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
3012 *
constraints_dirichlet_and_hanging_nodes_dual,
3020 *
constraints_dirichlet_and_hanging_nodes_dual.close();
3025 *
void DualSolver<dim>::assemble_matrix()
3032 *
const unsigned int dofs_per_cell_dual = fe_dual.dofs_per_cell;
3033 *
const unsigned int n_q_points = quadrature_formula.size();
3037 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell_dual);
3040 *
cell_dual = dof_handler_dual.begin_active(),
3041 *
endc_dual = dof_handler_dual.end(),
3042 *
cell = dof_handler.begin_active();
3046 *
for (; cell_dual != endc_dual; ++cell_dual, ++cell)
3047 *
if (cell_dual->is_locally_owned())
3049 *
fe_values.reinit(cell);
3051 *
fe_values_dual.reinit(cell_dual);
3054 *
std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3055 *
fe_values[displacement].get_function_symmetric_gradients(solution,
3058 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3062 *
constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3063 *
stress_strain_tensor_linearized,
3064 *
stress_strain_tensor);
3066 *
for (
unsigned int i = 0; i < dofs_per_cell_dual; ++i)
3069 *
stress_phi_i = stress_strain_tensor_linearized
3070 *
* fe_values_dual[displacement].symmetric_gradient(i, q_point);
3072 *
for (
unsigned int j = 0; j < dofs_per_cell_dual; ++j)
3074 *
* fe_values_dual[displacement].symmetric_gradient(j, q_point)
3075 *
* fe_values_dual.JxW(q_point));
3081 *
cell_dual->get_dof_indices(local_dof_indices);
3082 *
constraints_dirichlet_and_hanging_nodes_dual.distribute_local_to_global(cell_matrix,
3083 *
local_dof_indices,
3084 *
system_matrix_dual);
3092 *
void DualSolver<dim>::assemble_rhs()
3094 *
dual_functional->assemble_rhs (dof_handler, solution, constitutive_law,
3095 *
dof_handler_dual, system_rhs_dual);
3096 *
constraints_dirichlet_and_hanging_nodes_dual.condense (system_rhs_dual);
3101 *
void DualSolver<dim>::solve()
3105 * +++ direct solver +++++++++
3113 * After the decomposition, we can use A_direct like a
matrix representing
3114 * the inverse of our system
matrix, so to compute the solution we just
3115 * have to multiply with the right hand side vector:
3118 *
A_direct.vmult(solution_dual, system_rhs_dual);
3122 * ++++ iterative solver ++ CG ++++ doesn
't work
3123 * SolverControl solver_control (5000, 1e-12);
3124 * SolverCG<> cg (solver_control);
3128 * PreconditionSSOR<> preconditioner;
3129 * preconditioner.initialize(system_matrix_dual, 1.2);
3133 * cg.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3138 * ++++ iterative solver ++ BiCGStab ++++++ doesn't work
3145 * preconditioner.
initialize(system_matrix_dual, 1.0);
3149 * bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3154 * +++++++++++++++++++++++++++++++++++++++++++++++++
3160 *
constraints_dirichlet_and_hanging_nodes_dual.distribute (solution_dual);
3164 *
void DualSolver<dim>::output_results()
3166 *
std::string filename = (output_dir +
"dual-solution-" +
3168 *
std::ofstream output (filename.c_str());
3171 *
std::vector<std::string> solution_names;
3175 *
solution_names.push_back (
"displacement");
3178 *
solution_names.push_back (
"x_displacement");
3179 *
solution_names.push_back (
"y_displacement");
3182 *
solution_names.push_back (
"x_displacement");
3183 *
solution_names.push_back (
"y_displacement");
3184 *
solution_names.push_back (
"z_displacement");
3187 *
Assert (
false, ExcNotImplemented());
3189 *
data_out.add_data_vector (solution_dual, solution_names);
3190 *
data_out.build_patches ();
3191 *
data_out.write_vtk (output);
3195 *
void DualSolver<dim>::compute_error_DWR (
Vector<float> &estimated_error_per_cell)
3197 *
Assert (estimated_error_per_cell.size() == triangulation->n_global_active_cells(),
3198 *
ExcDimensionMismatch (estimated_error_per_cell.size(), triangulation->n_global_active_cells()));
3202 * solve the dual problem
3206 *
assemble_matrix ();
3209 *
output_results ();
3213 * compuate the dual weights
3220 *
constraints_dirichlet_and_hanging_nodes_dual,
3225 *
constraints_hanging_nodes);
3226 *
constraints_hanging_nodes.close();
3229 *
constraints_dirichlet_and_hanging_nodes_dual,
3232 *
constraints_hanging_nodes,
3237 * estimate the error
3247 *
const unsigned int n_q_points = quadrature_formula.size();
3248 *
std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3252 *
std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
3253 *
for (
unsigned int i=0; i!=n_q_points; ++i)
3255 *
cell_hessians[i].resize (dim);
3257 *
std::vector<Vector<double> > dual_weights_cell_values (n_q_points,
Vector<double>(dim));
3259 *
const EquationData::BodyForce<dim> body_force;
3260 *
std::vector<Vector<double> > body_force_values (n_q_points,
Vector<double>(dim));
3270 *
fe_face_values_neighbor (fe_dual, face_quadrature_formula,
3278 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
3279 *
std::vector<Vector<double> > jump_residual (n_face_q_points,
Vector<double>(dim));
3280 *
std::vector<Vector<double> > dual_weights_face_values (n_face_q_points,
Vector<double>(dim));
3282 *
std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
3283 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3285 *
cell_grads[i].resize (dim);
3287 *
std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
3288 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3290 *
neighbor_grads[i].resize (dim);
3298 *
typename std::map<typename DoFHandler<dim>::face_iterator,
Vector<double> >
3301 *
cell = dof_handler_dual.begin_active(),
3302 *
endc = dof_handler_dual.
end();
3303 *
for (; cell!=endc; ++cell)
3304 *
if (cell->is_locally_owned())
3306 *
for (
unsigned int face_no=0;
3307 *
face_no<GeometryInfo<dim>::faces_per_cell;
3310 *
face_integrals[cell->face(face_no)].reinit (dim);
3311 *
face_integrals[cell->face(face_no)] = -1e20;
3315 *
std::vector<Vector<float> > error_indicators_vector;
3316 *
error_indicators_vector.resize( triangulation->n_active_cells(),
3321 * ----------------- estimate_some -------------------------
3324 *
cell = dof_handler_dual.begin_active();
3325 *
unsigned int present_cell = 0;
3326 *
for (; cell!=endc; ++cell, ++present_cell)
3327 *
if (cell->is_locally_owned())
3331 * --------------- integrate_over_cell -------------------
3334 *
fe_values.reinit(cell);
3335 *
body_force.vector_value_list(fe_values.get_quadrature_points(),
3336 *
body_force_values);
3337 *
fe_values[displacement].get_function_symmetric_gradients(primal_solution,
3339 *
fe_values.get_function_hessians(primal_solution, cell_hessians);
3341 *
fe_values.get_function_values(dual_weights,
3342 *
dual_weights_cell_values);
3344 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3346 *
constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3347 *
stress_strain_tensor_linearized,
3348 *
stress_strain_tensor);
3349 *
constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
3350 *
cell_hessians[q_point],
3351 *
stress_strain_tensor_grad);
3353 *
for (
unsigned int i=0; i!=dim; ++i)
3355 *
error_indicators_vector[present_cell](i) +=
3356 *
body_force_values[q_point](i)*
3357 *
dual_weights_cell_values[q_point](i)*
3358 *
fe_values.JxW(q_point);
3359 *
for (
unsigned int j=0; j!=dim; ++j)
3361 *
for (
unsigned int k=0; k!=dim; ++k)
3363 *
for (
unsigned int l=0;
l!=dim; ++
l)
3365 *
error_indicators_vector[present_cell](i) +=
3366 *
( stress_strain_tensor[i][j][k][l]*
3367 *
0.5*(cell_hessians[q_point][k][l][j]
3369 *
cell_hessians[q_point][l][k][j])
3370 *
+ stress_strain_tensor_grad[i][j][k][
l][j] * strain_tensor[q_point][k][
l]
3372 *
dual_weights_cell_values[q_point](i) *
3373 *
fe_values.JxW(q_point);
3383 * -------------------------------------------------------
3384 * compute face_integrals
3387 *
for (
unsigned int face_no=0;
3388 *
face_no<GeometryInfo<dim>::faces_per_cell;
3391 *
if (cell->face(face_no)->at_boundary())
3393 *
for (
unsigned int id=0;
id!=dim; ++id)
3395 *
face_integrals[cell->face(face_no)](id) = 0;
3400 *
if ((cell->neighbor(face_no)->has_children() ==
false) &&
3401 *
(cell->neighbor(face_no)->level() == cell->level()) &&
3402 *
(cell->neighbor(face_no)->index() < cell->index()))
3405 *
if (cell->at_boundary(face_no) ==
false)
3406 *
if (cell->neighbor(face_no)->level() < cell->level())
3410 *
if (cell->face(face_no)->has_children() ==
false)
3414 * ------------- integrate_over_regular_face -----------
3417 *
fe_face_values_cell.reinit(cell, face_no);
3418 *
fe_face_values_cell.get_function_gradients (primal_solution,
3422 *
ExcInternalError());
3423 *
const unsigned int
3424 *
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3426 *
neighbor = cell->neighbor(face_no);
3428 *
fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
3429 *
fe_face_values_neighbor.get_function_gradients (primal_solution,
3432 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3434 *
q_cell_strain_tensor = 0.;
3435 *
q_neighbor_strain_tensor = 0.;
3436 *
for (
unsigned int i=0; i!=dim; ++i)
3438 *
for (
unsigned int j=0; j!=dim; ++j)
3440 *
q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3441 *
cell_grads[q_point][j][i] );
3442 *
q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3443 *
neighbor_grads[q_point][j][i] );
3447 *
constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3448 *
cell_stress_strain_tensor);
3449 *
constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3450 *
neighbor_stress_strain_tensor);
3452 *
jump_residual[q_point] = 0.;
3453 *
for (
unsigned int i=0; i!=dim; ++i)
3455 *
for (
unsigned int j=0; j!=dim; ++j)
3457 *
for (
unsigned int k=0; k!=dim; ++k)
3459 *
for (
unsigned int l=0;
l!=dim; ++
l)
3461 *
jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
3462 *
q_cell_strain_tensor[k][l]
3464 *
neighbor_stress_strain_tensor[i][j][k][l]*
3465 *
q_neighbor_strain_tensor[k][l] )*
3466 *
fe_face_values_cell.normal_vector(q_point)[j];
3474 *
fe_face_values_cell.get_function_values (dual_weights,
3475 *
dual_weights_face_values);
3478 *
face_integral_vector = 0;
3479 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3481 *
for (
unsigned int i=0; i!=dim; ++i)
3483 *
face_integral_vector(i) += jump_residual[q_point](i) *
3484 *
dual_weights_face_values[q_point](i) *
3485 *
fe_face_values_cell.JxW(q_point);
3489 *
Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
3490 *
ExcInternalError());
3492 *
for (
unsigned int i=0; i!=dim; ++i)
3494 *
Assert (face_integrals[cell->face(face_no)](i) == -1e20,
3495 *
ExcInternalError());
3496 *
face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
3502 * -----------------------------------------------------
3510 * ------------- integrate_over_irregular_face ---------
3514 *
face = cell->face(face_no);
3516 *
neighbor = cell->neighbor(face_no);
3518 *
ExcInternalError());
3519 *
Assert (neighbor->has_children(),
3520 *
ExcInternalError());
3522 *
const unsigned int
3523 *
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3525 *
for (
unsigned int subface_no=0;
3526 *
subface_no<face->n_children(); ++subface_no)
3529 *
neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
3530 *
Assert (neighbor_child->face(neighbor_neighbor) ==
3531 *
cell->face(face_no)->child(subface_no),
3532 *
ExcInternalError());
3534 *
fe_subface_values_cell.reinit (cell, face_no, subface_no);
3535 *
fe_subface_values_cell.get_function_gradients (primal_solution,
3537 *
fe_face_values_neighbor.reinit (neighbor_child,
3538 *
neighbor_neighbor);
3539 *
fe_face_values_neighbor.get_function_gradients (primal_solution,
3542 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3544 *
q_cell_strain_tensor = 0.;
3545 *
q_neighbor_strain_tensor = 0.;
3546 *
for (
unsigned int i=0; i!=dim; ++i)
3548 *
for (
unsigned int j=0; j!=dim; ++j)
3550 *
q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3551 *
cell_grads[q_point][j][i] );
3552 *
q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3553 *
neighbor_grads[q_point][j][i] );
3557 *
constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3558 *
cell_stress_strain_tensor);
3559 *
constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3560 *
neighbor_stress_strain_tensor);
3562 *
jump_residual[q_point] = 0.;
3563 *
for (
unsigned int i=0; i!=dim; ++i)
3565 *
for (
unsigned int j=0; j!=dim; ++j)
3567 *
for (
unsigned int k=0; k!=dim; ++k)
3569 *
for (
unsigned int l=0;
l!=dim; ++
l)
3571 *
jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
3572 *
q_cell_strain_tensor[k][l]
3574 *
neighbor_stress_strain_tensor[i][j][k][l]*
3575 *
q_neighbor_strain_tensor[k][l] )*
3576 *
fe_face_values_neighbor.normal_vector(q_point)[j];
3584 *
fe_face_values_neighbor.get_function_values (dual_weights,
3585 *
dual_weights_face_values);
3588 *
face_integral_vector = 0;
3589 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3591 *
for (
unsigned int i=0; i!=dim; ++i)
3593 *
face_integral_vector(i) += jump_residual[q_point](i) *
3594 *
dual_weights_face_values[q_point](i) *
3595 *
fe_face_values_neighbor.JxW(q_point);
3599 *
for (
unsigned int i=0; i!=dim; ++i)
3601 *
face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
3608 *
for (
unsigned int subface_no=0;
3609 *
subface_no<face->n_children(); ++subface_no)
3611 *
Assert (face_integrals.find(face->child(subface_no)) !=
3612 *
face_integrals.end(),
3613 *
ExcInternalError());
3614 *
for (
unsigned int i=0; i!=dim; ++i)
3616 *
Assert (face_integrals[face->child(subface_no)](i) != -1e20,
3617 *
ExcInternalError());
3618 *
sum(i) += face_integrals[face->child(subface_no)](i);
3621 *
for (
unsigned int i=0; i!=dim; ++i)
3623 *
face_integrals[face](i) =
sum(i);
3629 * -----------------------------------------------------
3639 * ----------------------------------------------------------
3646 *
cell = dof_handler_dual.begin_active();
3647 *
for (; cell!=endc; ++cell, ++present_cell)
3648 *
if (cell->is_locally_owned())
3650 *
for (
unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
3653 *
Assert(face_integrals.find(cell->face(face_no)) !=
3654 *
face_integrals.end(),
3655 *
ExcInternalError());
3657 *
for (
unsigned int id=0;
id!=dim; ++id)
3659 *
error_indicators_vector[present_cell](id)
3660 *
-= 0.5*face_integrals[cell->face(face_no)](id);
3665 *
estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
3675 * <a name=
"elastoplastic.cc-ThecodePlasticityContactProblemcodeclasstemplate"></a>
3676 * <h3>The <code>PlasticityContactProblem</code>
class template</h3>
3680 * This is the
main class of this program and supplies all
functions
3681 * and variables needed to describe
3682 * the nonlinear contact problem. It is
3683 * close to @ref step_41
"step-41" but with some additional
3684 * features like handling hanging nodes,
3685 * a Newton method,
using Trilinos and p4est
3686 *
for parallel distributed computing.
3687 * To deal with hanging nodes makes
3688 * life a bit more complicated since
3690 * We create a Newton method
for the
3691 * active
set method
for the contact
3692 * situation and to handle the nonlinear
3693 *
operator for the constitutive law.
3697 * The
general layout of
this class is very much like
for most other tutorial programs.
3698 * To make our life a bit easier,
this class reads a
set of input parameters from an input file. These
3699 * parameters,
using the
ParameterHandler class, are declared in the <code>declare_parameters</code>
3700 * function (which is
static so that it can be called before we even create an
object of the current
3701 * type), and a
ParameterHandler object that has been used to read an input file will then be passed
3702 * to the constructor of
this class.
3706 * The remaining member
functions are by and large as we have seen in several of the other tutorial
3707 * programs, though with additions
for the current nonlinear system. We will comment on their purpose
3708 * as we get to them further below.
3711 *
template <
int dim>
3712 *
class ElastoPlasticProblem
3722 *
void make_grid ();
3723 *
void setup_system ();
3724 *
void compute_dirichlet_constraints ();
3728 *
void solve_newton_system ();
3729 *
void solve_newton ();
3730 *
void compute_error ();
3732 *
void refine_grid ();
3734 *
void output_results (
const std::string &filename_base);
3738 * Next are three
functions that handle the history variables stored in each
3740 *
set up a pristine state
for the history variables. It only works on
3741 * those quadrature points on cells that belong to the present processor:
3744 *
void setup_quadrature_point_history ();
3748 * The
second one updates the history variables at the
end of each
3752 *
void update_quadrature_point_history ();
3756 * As far as member variables are concerned, we start with ones that we use to
3757 * indicate the
MPI universe
this program runs on, and then two
numbers
3758 * telling us how many participating processors there are, and where in
3759 *
this world we are., a stream we use to let
3760 * exactly
one processor produce output to the console (see @ref step_17
"step-17") and
3761 * a variable that is used to time the various sections of the program:
3772 * The next
group describes the mesh and the finite element space.
3773 * In particular,
for this parallel program, the finite element
3774 * space has associated with it variables that indicate which degrees
3775 * of freedom live on the current processor (the index sets, see
3776 * also @ref step_40
"step-40" and the @ref distributed documentation module) as
3777 * well as a variety of constraints: those imposed by hanging nodes,
3778 * by Dirichlet boundary conditions, and by the active
set of
3780 * here, the
first only contains hanging node constraints, the
3781 *
second also those associated with Dirichlet boundary conditions,
3782 * and the third these plus the contact constraints.
3786 * The variable <code>active_set</code> consists of those degrees
3787 * of freedom constrained by the contact, and we use
3788 * <code>fraction_of_plastic_q_points_per_cell</code> to keep
3789 * track of the fraction of quadrature points on each cell where
3790 * the stress equals the yield stress. The latter is only used to
3791 * create graphical output showing the plastic zone, but not
for
3792 * any further computation; the variable is a member variable of
3793 *
this class since the information is computed as a by-product
3794 * of computing the residual, but is used only much later. (Note
3795 * that the vector is a vector of length
equal to the number of
3796 * active cells on the <i>local mesh</i>; it is never used to
3797 * exchange information between processors and can therefore be
3798 * a regular deal.II vector.)
3801 *
const unsigned int n_initial_global_refinements;
3804 *
const unsigned int fe_degree;
3818 * One difference of
this program is that we declare the quadrature
3819 * formula in the
class declaration. The reason is that in all the other
3820 * programs, it didn
't do much harm if we had used different quadrature
3821 * formulas when computing the matrix and the right hand side, for
3822 * example. However, in the present case it does: we store information in
3823 * the quadrature points, so we have to make sure all parts of the program
3824 * agree on where they are and how many there are on each cell. Thus, let
3825 * us first declare the quadrature formula that will be used throughout...
3828 * const QGauss<dim> quadrature_formula;
3829 * const QGauss<dim - 1> face_quadrature_formula;
3833 * ... and then also have a vector of history objects, one per quadrature
3834 * point on those cells for which we are responsible (i.e. we don't store
3835 * history
data for quadrature points on cells that are owned by other
3839 *
std::vector<PointHistory<dim> > quadrature_point_history;
3843 * The way
this object is accessed is through a <code>user
pointer</code>
3844 * that each cell, face, or edge holds: it is a <code>
void*</code>
pointer
3845 * that can be used by application programs to associate arbitrary
data to
3846 * cells, faces, or edges. What the program actually does with
this data
3847 * is within its own responsibility, the library just allocates some space
3848 *
for these pointers, and application programs can
set and read the
3849 * pointers
for each of these objects.
3856 * The next block of variables corresponds to the solution
3857 * and the linear systems we need to form. In particular,
this
3858 * includes the Newton
matrix and right hand side; the vector
3859 * that corresponds to the residual (i.e., the Newton right hand
3860 * side) but from which we have not eliminated the various
3861 * constraints and that is used to determine which degrees of
3862 * freedom need to be constrained in the next iteration; and
3863 * a vector that corresponds to the
diagonal of the @f$B@f$
matrix
3864 * briefly mentioned in the introduction and discussed in the
3865 * accompanying paper.
3877 * The next block of variables is then related to the time dependent
3878 * nature of the problem: they denote the length of the time interval
3879 * which we want to simulate, the present time and number of time step,
3880 * and length of present timestep:
3883 *
double present_time;
3884 *
double present_timestep;
3886 *
unsigned int timestep_no;
3890 * The next block contains the variables that describe the material
3894 *
const double e_modulus, nu, sigma_0,
gamma;
3895 *
ConstitutiveLaw<dim> constitutive_law;
3899 * And then there is an assortment of other variables that are used
3900 * to identify the mesh we are asked to build as selected by the
3901 * parameter file, the obstacle that is being pushed into the
3902 * deformable body, the mesh refinement strategy, whether to transfer
3903 * the solution from
one mesh to the next, and how many mesh
3904 * refinement cycles to perform. As possible, we mark these kinds
3905 * of variables as <code>
const</code> to help the reader identify
3906 * which ones may or may not be modified later on (the output directory
3907 * being an exception -- it is never modified outside the constructor
3908 * but it is awkward to initialize in the member-initializer-list
3909 * following the colon in the constructor since there we have only
3910 * one shot at setting it; the same is
true for the mesh refinement
3916 *
struct RefinementStrategy
3921 *
refine_percentage,
3925 *
typename RefinementStrategy::value refinement_strategy;
3927 *
struct ErrorEstimationStrategy
3933 *
weighted_residual_error,
3934 *
weighted_kelly_error
3937 *
typename ErrorEstimationStrategy::value error_estimation_strategy;
3941 *
const bool transfer_solution;
3942 *
std::string output_dir;
3947 *
unsigned int current_refinement_cycle;
3949 *
const double max_relative_error;
3950 *
float relative_error;
3952 *
const bool show_stresses;
3959 * <a name=
"elastoplastic.cc-ImplementationofthecodePlasticityContactProblemcodeclass"></a>
3960 * <h3>Implementation of the <code>PlasticityContactProblem</code>
class</h3>
3965 * <a name=
"elastoplastic.cc-PlasticityContactProblemdeclare_parameters"></a>
3966 * <h4>PlasticityContactProblem::declare_parameters</h4>
3970 * Let us start with the declaration of
run-time parameters that can be
3971 * selected in the input file. These
values will be read back in the
3972 * constructor of
this class to initialize the member variables of this
3976 *
template <
int dim>
3980 *
prm.declare_entry(
"polynomial degree",
"1",
3982 *
"Polynomial degree of the FE_Q finite element space, typically 1 or 2.");
3983 *
prm.declare_entry(
"number of initial refinements",
"2",
3985 *
"Number of initial global mesh refinement steps before "
3986 *
"the first computation.");
3987 *
prm.declare_entry(
"refinement strategy",
"percentage",
3989 *
"Mesh refinement strategy:\n"
3990 *
" global: one global refinement\n"
3991 *
" percentage: a fixed percentage of cells gets refined using the selected error estimator.");
3992 *
prm.declare_entry(
"error estimation strategy",
"kelly_error",
3994 *
"Error estimation strategy:\n"
3995 *
" kelly_error: Kelly error estimator\n"
3996 *
" residual_error: residual-based error estimator\n"
3997 *
" weighted_residual_error: dual weighted residual (Goal-oriented) error estimator.\n");
3998 *
prm.declare_entry(
"maximum relative error",
"0.05",
4000 *
"maximum relative error which plays the role of a criteria for refinement.");
4001 *
prm.declare_entry(
"number of cycles",
"5",
4003 *
"Number of adaptive mesh refinement cycles to run.");
4004 *
prm.declare_entry(
"output directory",
"",
4006 *
"Directory for output files (graphical output and benchmark "
4007 *
"statistics). If empty, use the current directory.");
4008 *
prm.declare_entry(
"transfer solution",
"true",
4010 *
"Whether the solution should be used as a starting guess "
4011 *
"for the next finer mesh. If false, then the iteration starts at "
4012 *
"zero on every mesh.");
4013 *
prm.declare_entry(
"base mesh",
"Thick_tube_internal_pressure",
4015 *
"Perforated_strip_tension|Cantiliver_beam_3d"),
4016 *
"Select the shape of the domain: 'box' or 'half sphere'");
4017 *
prm.declare_entry(
"elasticity modulus",
"2.e11",
4019 *
"Elasticity modulus of the material in MPa (N/mm2)");
4020 *
prm.declare_entry(
"Poissons ratio",
"0.3",
4022 *
"Poisson's ratio of the material");
4023 *
prm.declare_entry(
"yield stress",
"2.e11",
4025 *
"Yield stress of the material in MPa (N/mm2)");
4026 *
prm.declare_entry(
"isotropic hardening parameter",
"0.",
4028 *
"Isotropic hardening parameter of the material");
4029 *
prm.declare_entry(
"show stresses",
"false",
4031 *
"Whether illustrates the stresses and von Mises stresses or not.");
4040 * <a name=
"elastoplastic.cc-ThecodePlasticityContactProblemcodeconstructor"></a>
4041 * <h4>The <code>PlasticityContactProblem</code> constructor</h4>
4045 * Given the declarations of member variables as well as the
4046 * declarations of
run-time parameters that are read from the input
4047 * file, there is
nothing surprising in
this constructor. In the body
4048 * we initialize the mesh refinement strategy and the output directory,
4049 * creating such a directory
if necessary.
4052 *
template <
int dim>
4053 *
ElastoPlasticProblem<dim>::
4056 *
mpi_communicator(MPI_COMM_WORLD),
4059 *
pcout(std::cout, this_mpi_process == 0),
4063 *
n_initial_global_refinements (prm.get_integer(
"number of initial refinements")),
4064 *
triangulation(mpi_communicator),
4065 *
fe_degree (prm.get_integer(
"polynomial degree")),
4067 *
dof_handler(triangulation),
4068 *
quadrature_formula (fe_degree + 1),
4069 *
face_quadrature_formula (fe_degree + 1),
4071 *
e_modulus (prm.get_double(
"elasticity modulus")),
4072 *
nu (prm.get_double(
"Poissons ratio")),
4073 *
sigma_0(prm.get_double(
"yield stress")),
4074 *
gamma (prm.get_double(
"isotropic hardening parameter")),
4075 *
constitutive_law (e_modulus,
4080 *
base_mesh (prm.get(
"base mesh")),
4082 *
transfer_solution (prm.get_bool(
"transfer solution")),
4084 *
table_results_2(),
4085 *
table_results_3(),
4086 *
max_relative_error (prm.get_double(
"maximum relative error")),
4087 *
show_stresses (prm.get_bool(
"show stresses"))
4089 *
std::string strat = prm.get(
"refinement strategy");
4090 *
if (strat ==
"global")
4091 *
refinement_strategy = RefinementStrategy::refine_global;
4092 *
else if (strat ==
"percentage")
4093 *
refinement_strategy = RefinementStrategy::refine_percentage;
4097 *
strat = prm.get(
"error estimation strategy");
4098 *
if (strat ==
"kelly_error")
4099 *
error_estimation_strategy = ErrorEstimationStrategy::kelly_error;
4100 *
else if (strat ==
"residual_error")
4101 *
error_estimation_strategy = ErrorEstimationStrategy::residual_error;
4102 *
else if (strat ==
"weighted_residual_error")
4103 *
error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error;
4107 *
output_dir = prm.get(
"output directory");
4108 *
if (output_dir !=
"" && *(output_dir.rbegin()) !=
'/')
4109 *
output_dir +=
"/";
4110 *
mkdir(output_dir.c_str(), 0777);
4112 *
pcout <<
" Using output directory '" << output_dir <<
"'" << std::endl;
4113 *
pcout <<
" FE degree " << fe_degree << std::endl;
4114 *
pcout <<
" transfer solution "
4115 *
<< (transfer_solution ?
"true" :
"false") << std::endl;
4123 * <a name=
"elastoplastic.cc-PlasticityContactProblemmake_grid"></a>
4124 * <h4>PlasticityContactProblem::make_grid</h4>
4128 * The next block deals with constructing the starting mesh.
4129 * We will use the following helper function and the
first
4130 * block of the <code>make_grid()</code> to construct a
4131 * mesh that corresponds to a half sphere. deal.II has a function
4132 * that creates such a mesh, but it is in the wrong location
4133 * and facing the wrong direction, so we need to
shift and
rotate
4134 * it a bit before
using it.
4138 * For later
reference, as described in the documentation of
4140 * has boundary indicator
zero,
while the remainder has boundary
4145 *
rotate_half_sphere (
const Point<3> &in)
4147 *
return Point<3>(in(2), in(1), -in(0));
4150 *
template <
int dim>
4152 *
ElastoPlasticProblem<dim>::make_grid ()
4154 *
if (base_mesh ==
"Timoshenko beam")
4158 *
const double length = .48,
4162 *
point_2(length, depth/2);
4164 *
std::vector<unsigned int> repetitions(2);
4165 *
repetitions[0] = 4;
4166 *
repetitions[1] = 1;
4172 * give the indicators to boundaries
for specification,
4179 * |________________|
4181 * 0 to essential boundary conditions (left edge) which are as
default
4182 * 100 to the null boundaries (upper and lower edges) where we
do not need to take care of them
4183 * 5 to the natural boundaries (right edge)
for imposing the traction force
4187 *
cell = triangulation.begin_active(),
4188 *
endc = triangulation.end();
4189 *
for (; cell!=endc; ++cell)
4193 *
if ( std::fabs(cell->face(face)->center()(0)-length) < 1e-12 )
4195 *
cell->face(face)->set_manifold_id(5);
4197 *
else if ( ( std::fabs(cell->face(face)->center()(1)-(depth/2)) < 1e-12 )
4199 *
( std::fabs(cell->face(face)->center()(1)-(-depth/2)) < 1e-12 ) )
4201 *
cell->face(face)->set_manifold_id(100);
4207 *
triangulation.refine_global(n_initial_global_refinements);
4210 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4214 * Example 1 from the paper: Zhong Z., ....
A new numerical method
for determining
4215 * collapse
load-carrying capacity of structure made of elasto-plastic material,
4216 *
J. Cent. South Univ. (2014) 21: 398-404
4222 *
const double inner_radius = .1,
4230 * give the indicators to boundaries
for specification,
4248 * 0 - inner boundary - natural boundary condition - impose the traction force
4249 * 1 - outer boundary - free boundary - we
do not need to take care of them
4250 * 2 - left boundary - essential boundary condition - constrained to move along the x direction
4251 * 3 - bottom boundary - essential boundary condition - constrained to move along the y direction
4258 *
triangulation.set_manifold (0, inner_boundary_description);
4261 *
triangulation.set_manifold (1, outer_boundary_description);
4263 *
triangulation.refine_global(n_initial_global_refinements);
4265 *
triangulation.reset_manifold (0);
4266 *
triangulation.reset_manifold (1);
4269 *
else if (base_mesh ==
"Perforated_strip_tension")
4273 * Example 2 from the paper: Zhong Z., ....
A new numerical method
for determining
4274 * collapse
load-carrying capacity of structure made of elasto-plastic material,
4275 *
J. Cent. South Univ. (2014) 21: 398-404
4280 *
const int dim_2d = 2;
4282 *
const double inner_radius = 0.05,
4285 *
thickness = 0.004;
4298 *
const double eps = 1
e-7 * inner_radius;
4308 * Modify the triangulation_1
4312 *
cell = triangulation_1.begin_active(),
4313 *
endc = triangulation_1.end();
4314 *
std::vector<bool> treated_vertices(triangulation_1.n_vertices(),
false);
4315 *
for (; cell != endc; ++cell)
4317 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4318 *
if (cell->face(f)->at_boundary() && cell->face(f)->center()(0)>eps &&
4319 *
cell->face(f)->center()(1)>eps )
4323 * distance of the face center from the center
4326 *
point(0) = cell->face(f)->center()(0) - center_2d(0);
4327 *
point(1) = cell->face(f)->center()(1) - center_2d(1);
4330 *
for (
unsigned int v=0; v < GeometryInfo<dim_2d>::vertices_per_face; ++v)
4332 *
unsigned int vv = cell->face(f)->vertex_index(v);
4333 *
if (treated_vertices[vv] ==
false)
4335 *
treated_vertices[vv] =
true;
4351 * Make the triangulation_2, a rectangular above the triangulation_1
4364 * make the triangulation_2d and
refine it
4370 * Merge the two triangulation_1 and triangulation_2
4377 * Assign boundary indicators to the boundary faces
4398 *
cell = triangulation_2d.begin_active(),
4399 *
endc = triangulation_2d.end();
4400 *
for (; cell != endc; ++cell)
4402 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4404 *
if (cell->face(f)->at_boundary())
4406 *
if ( std::fabs(cell->face(f)->center()(1)) < eps )
4408 *
cell->face(f)->set_manifold_id(1);
4410 *
else if ( std::fabs(cell->face(f)->center()(0)-
outer_radius) < eps )
4412 *
cell->face(f)->set_manifold_id(2);
4414 *
else if ( std::fabs(cell->face(f)->center()(1)-
height) < eps )
4416 *
cell->face(f)->set_manifold_id(3);
4418 *
else if ( std::fabs(cell->face(f)->center()(0)) < eps )
4420 *
cell->face(f)->set_manifold_id(4);
4424 *
cell->face(f)->set_all_boundary_ids(10);
4434 *
triangulation_2d.set_manifold (10, inner_boundary_description);
4436 *
triangulation_2d.refine_global(3);
4438 *
triangulation_2d.reset_manifold (10);
4443 * Extrude the triangulation_2d and make it 3
d
4445 * 2, thickness, triangulation);
4449 *
2, thickness, triangulation);
4453 * Assign boundary indicators to the boundary faces
4474 *
Point<dim> center(center_2d(0), center_2d(1), 0);
4477 *
cell = triangulation.begin_active(),
4478 *
endc = triangulation.end();
4479 *
for (; cell != endc; ++cell)
4481 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4483 *
if (cell->face(f)->at_boundary())
4485 *
dist_vector = cell->face(f)->center() - center;
4487 *
if ( std::fabs(dist_vector[1]) < eps )
4489 *
cell->face(f)->set_manifold_id(1);
4491 *
else if ( std::fabs(dist_vector[0]-
outer_radius) < eps )
4493 *
cell->face(f)->set_manifold_id(2);
4495 *
else if ( std::fabs(dist_vector[1]-
height) < eps )
4497 *
cell->face(f)->set_manifold_id(3);
4499 *
else if ( std::fabs(dist_vector[0]) < eps )
4501 *
cell->face(f)->set_manifold_id(4);
4503 *
else if ( std::fabs(dist_vector[2]) < eps )
4505 *
cell->face(f)->set_manifold_id(5);
4507 *
else if ( std::fabs(dist_vector[2]-thickness) < eps )
4509 *
cell->face(f)->set_manifold_id(6);
4513 *
cell->face(f)->set_all_boundary_ids(10);
4523 *
triangulation.set_manifold (10, inner_boundary_description);
4525 *
triangulation.refine_global(n_initial_global_refinements);
4527 *
triangulation.reset_manifold (10);
4530 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4534 *
A rectangular tube made of Aluminium
4536 * approximation of beam 17250
4537 * units are in meter
4545 *
const int dim_2d = 2;
4547 *
const double length = .7,
4550 *
thickness_web = 10
e-3,
4551 *
thickness_flange = 10
e-3;
4559 *
const double eps = 1
e-7 * width;
4562 * Make the triangulation_b, a rectangular at the bottom of rectangular tube
4567 *
point2 (width/2, -(
height/2)+thickness_flange);
4569 *
std::vector<unsigned int> repetitions(dim_2d);
4570 *
repetitions[0] = 8;
4571 *
repetitions[1] = 1;
4578 * Make the triangulation_t, a rectangular at the top of rectangular tube
4583 *
point2 (width/2,
height/2);
4585 *
std::vector<unsigned int> repetitions(dim_2d);
4586 *
repetitions[0] = 8;
4587 *
repetitions[1] = 1;
4594 * Make the triangulation_l, a rectangular at the left of rectangular tube
4599 *
point2 (-(width/2)+thickness_web, (
height/2)-thickness_flange);
4601 *
std::vector<unsigned int> repetitions(dim_2d);
4602 *
repetitions[0] = 1;
4603 *
repetitions[1] = 18;
4610 * Make the triangulation_r, a rectangular at the right of rectangular tube
4615 *
point2 (width/2, (
height/2)-thickness_flange);
4617 *
std::vector<unsigned int> repetitions(dim_2d);
4618 *
repetitions[0] = 1;
4619 *
repetitions[1] = 18;
4626 * make the triangulation_2d
4632 * merging every two triangles to make triangulation_2d
4636 *
triangulation_blr;
4645 * Extrude the triangulation_2d and make it 3
d
4648 *
const unsigned int n_slices =
static_cast<int>(length*1000/20) + 1;
4650 *
n_slices, length, triangulation);
4654 * Assign boundary indicators to the boundary faces
4696 *
cell = triangulation.begin_active(),
4697 *
endc = triangulation.end();
4698 *
for (; cell != endc; ++cell)
4700 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4702 *
if (cell->face(f)->at_boundary())
4704 *
dist_vector = cell->face(f)->center() - center;
4706 *
if ( std::fabs(dist_vector[2]) < eps )
4708 *
cell->face(f)->set_manifold_id(1);
4710 *
else if ( std::fabs(dist_vector[1]-(
height/2)) < eps )
4712 *
cell->face(f)->set_manifold_id(2);
4716 *
cell->face(f)->set_all_boundary_ids(0);
4725 *
triangulation.refine_global(n_initial_global_refinements);
4733 *
pcout <<
" Number of active cells: "
4734 *
<< triangulation.n_active_cells()
4743 * <a name=
"elastoplastic.cc-PlasticityContactProblemsetup_system"></a>
4744 * <h4>PlasticityContactProblem::setup_system</h4>
4748 * The next piece in the puzzle is to
set up the
DoFHandler, resize
4749 * vectors and take care of various other status variables such as
4750 *
index sets and constraint matrices.
4754 * In the following, each
group of operations is put into a brace-enclosed
4755 * block that is being timed by the variable declared at the top of the
4757 * timed section, the destructor that is called at the
end of the block
4761 *
template <
int dim>
4763 *
ElastoPlasticProblem<dim>::setup_system ()
4769 *
dof_handler.distribute_dofs(fe);
4770 *
pcout <<
" Number of degrees of freedom: "
4771 *
<< dof_handler.n_dofs()
4774 *
locally_owned_dofs = dof_handler.locally_owned_dofs();
4775 *
locally_relevant_dofs =
4782 *
constraints_hanging_nodes.reinit(locally_relevant_dofs);
4784 *
constraints_hanging_nodes);
4785 *
constraints_hanging_nodes.close();
4787 *
pcout <<
" Number of active cells: "
4788 *
<< triangulation.n_global_active_cells() << std::endl
4789 *
<<
" Number of degrees of freedom: " << dof_handler.n_dofs()
4792 *
compute_dirichlet_constraints();
4798 *
if (timestep_no==1 || current_refinement_cycle!=0)
4800 *
solution.reinit(locally_relevant_dofs, mpi_communicator);
4802 *
incremental_displacement.reinit(locally_relevant_dofs, mpi_communicator);
4803 *
newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
4804 *
newton_rhs_residual.reinit(locally_owned_dofs, mpi_communicator);
4805 *
fraction_of_plastic_q_points_per_cell.reinit(triangulation.n_active_cells());
4810 * Finally, we
set up sparsity patterns and matrices.
4811 * We temporarily (ab)use the system
matrix to also build the (diagonal)
4812 *
matrix that we use in eliminating degrees of freedom that are in contact
4813 * with the obstacle, but we then immediately
set the Newton
matrix back
4820 *
mpi_communicator);
4823 *
constraints_dirichlet_and_hanging_nodes,
false,
4824 *
this_mpi_process);
4826 *
newton_matrix.reinit(sp);
4834 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_dirichlet_constraints"></a>
4835 * <h4>PlasticityContactProblem::compute_dirichlet_constraints</h4>
4839 * This function, broken out of the preceding
one, computes the constraints
4840 * associated with Dirichlet-type boundary conditions and puts them into the
4841 * <code>constraints_dirichlet_and_hanging_nodes</code> variable by merging
4842 * with the constraints that come from hanging nodes.
4846 * As laid out in the introduction, we need to distinguish between two
4848 * - If the domain is a box, we
set the displacement to
zero at the bottom,
4849 * and allow vertical movement in z-direction along the sides. As
4850 * shown in the <code>make_grid()</code> function, the former corresponds
4851 * to boundary indicator 6, the latter to 8.
4852 * - If the domain is a half sphere, then we impose
zero displacement along
4853 * the curved part of the boundary, associated with boundary indicator
zero.
4856 *
template <
int dim>
4858 *
ElastoPlasticProblem<dim>::compute_dirichlet_constraints ()
4860 *
constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
4861 *
constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
4863 *
std::vector<bool> component_mask(dim);
4865 *
if (base_mesh ==
"Timoshenko beam")
4869 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4870 *
constraints_dirichlet_and_hanging_nodes,
4873 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4877 * the boundary x = 0
4880 *
component_mask[0] =
true;
4881 *
component_mask[1] =
false;
4884 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4885 *
constraints_dirichlet_and_hanging_nodes,
4889 * the boundary y = 0
4892 *
component_mask[0] =
false;
4893 *
component_mask[1] =
true;
4896 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4897 *
constraints_dirichlet_and_hanging_nodes,
4900 *
else if (base_mesh ==
"Perforated_strip_tension")
4904 * the boundary x = 0
4907 *
component_mask[0] =
true;
4908 *
component_mask[1] =
false;
4909 *
component_mask[2] =
false;
4912 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4913 *
constraints_dirichlet_and_hanging_nodes,
4917 * the boundary y = 0
4920 *
component_mask[0] =
false;
4921 *
component_mask[1] =
true;
4922 *
component_mask[2] =
false;
4925 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4926 *
constraints_dirichlet_and_hanging_nodes,
4930 * the boundary y = imposed incremental displacement
4933 *
component_mask[0] =
false;
4934 *
component_mask[1] =
true;
4935 *
component_mask[2] =
false;
4938 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4939 *
constraints_dirichlet_and_hanging_nodes,
4942 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4946 * the boundary x = y = z = 0
4949 *
component_mask[0] =
true;
4950 *
component_mask[1] =
true;
4951 *
component_mask[2] =
true;
4954 *
EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4955 *
constraints_dirichlet_and_hanging_nodes,
4964 *
constraints_dirichlet_and_hanging_nodes.close();
4971 * <a name=
"elastoplastic.cc-PlasticityContactProblemassemble_newton_system"></a>
4972 * <h4>PlasticityContactProblem::assemble_newton_system</h4>
4976 * Given the complexity of the problem, it may come as a bit of a surprise
4977 * that assembling the linear system we have to solve in each Newton iteration
4978 * is actually fairly straightforward. The following function builds the Newton
4979 * right hand side and Newton
matrix. It looks fairly innocent because the
4980 * heavy lifting happens in the call to
4981 * <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
4983 * constraints we have previously computed.
4986 *
template <
int dim>
4988 *
ElastoPlasticProblem<dim>::
4995 *
if (base_mesh ==
"Timoshenko beam")
4997 *
traction_surface_id = 5;
4999 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5001 *
traction_surface_id = 0;
5003 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5005 *
traction_surface_id = 2;
5019 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5020 *
const unsigned int n_q_points = quadrature_formula.size();
5021 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5024 *
const EquationData::BodyForce<dim> body_force;
5025 *
std::vector<Vector<double> > body_force_values(n_q_points,
5028 *
const EquationData::
5029 *
IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5030 *
std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5036 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5040 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5043 *
std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
5046 *
cell = dof_handler.begin_active(),
5047 *
endc = dof_handler.end();
5051 *
for (; cell != endc; ++cell)
5052 *
if (cell->is_locally_owned())
5054 *
fe_values.reinit(cell);
5058 *
fe_values[displacement].get_function_symmetric_gradients(delta_linearization_point,
5059 *
incremental_strain_tensor);
5063 * For assembling the local right hand side contributions, we need
5064 * to access the prior linearized stress
value in
this quadrature
5065 *
point. To get it, we use the user
pointer of
this cell that
5066 * points into the global array to the quadrature
point data
5067 * corresponding to the
first quadrature
point of the present cell,
5068 * and then add an offset corresponding to the
index of the
5069 * quadrature
point we presently consider:
5072 *
const PointHistory<dim> *local_quadrature_points_history
5073 *
=
reinterpret_cast<PointHistory<dim>*
>(cell->user_pointer());
5074 *
Assert (local_quadrature_points_history >=
5075 *
&quadrature_point_history.front(),
5076 *
ExcInternalError());
5077 *
Assert (local_quadrature_points_history <
5078 *
&quadrature_point_history.back(),
5079 *
ExcInternalError());
5083 * In addition, we need the
values of the external body forces at
5084 * the quadrature points on
this cell:
5087 *
body_force.vector_value_list(fe_values.get_quadrature_points(),
5088 *
body_force_values);
5090 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5093 *
tmp_strain_tensor_qpoint = local_quadrature_points_history[q_point].old_strain
5094 *
+ incremental_strain_tensor[q_point];
5098 *
constitutive_law.get_linearized_stress_strain_tensors(tmp_strain_tensor_qpoint,
5099 *
stress_strain_tensor_linearized,
5100 *
stress_strain_tensor);
5103 *
for (
unsigned int i = 0; i < dim; ++i)
5105 *
rhs_values_body_force[i] = body_force_values[q_point][i];
5108 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5112 * Having computed the stress-strain tensor and its linearization,
5113 * we can now put together the parts of the
matrix and right hand side.
5114 * In both, we need the linearized stress-strain tensor times the
5115 *
symmetric gradient of @f$\varphi_i@f$, i.e. the term @f$I_\Pi\varepsilon(\varphi_i)@f$,
5116 * so we introduce an abbreviation of
this term. Recall that the
5117 *
matrix corresponds to the bilinear form
5118 * @f$A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))@f$ in the
5119 * notation of the accompanying publication, whereas the right
5120 * hand side is @f$F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))@f$
5121 * where @f$u@f$ is the current linearization points (typically the last solution).
5122 * This might suggest that the right hand side will be
zero if the material
5123 * is completely elastic (where @f$I_\Pi=P_\Pi@f$) but
this ignores the fact
5124 * that the right hand side will also contain contributions from
5125 * non-homogeneous constraints due to the contact.
5129 * The code block that follows
this adds contributions that are due to
5130 * boundary forces, should there be any.
5134 *
stress_phi_i = stress_strain_tensor_linearized
5135 *
* fe_values[displacement].symmetric_gradient(i, q_point);
5137 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
5139 *
* fe_values[displacement].symmetric_gradient(j, q_point)
5140 *
* fe_values.JxW(q_point));
5144 *
* incremental_strain_tensor[q_point] )
5146 *
( ( stress_strain_tensor
5147 *
* fe_values[displacement].symmetric_gradient(i, q_point))
5148 *
* tmp_strain_tensor_qpoint )
5150 *
( fe_values[displacement].
value(i, q_point)
5151 *
* rhs_values_body_force )
5152 *
) * fe_values.JxW(q_point);
5157 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
5158 *
if (cell->face(face)->at_boundary()
5160 *
cell->face(face)->boundary_id() == traction_surface_id)
5162 *
fe_values_face.reinit(cell, face);
5164 *
boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5165 *
boundary_force_values);
5167 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
5170 *
for (
unsigned int i = 0; i < dim; ++i)
5172 *
rhs_values[i] = boundary_force_values[q_point][i];
5174 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5175 *
cell_rhs(i) += (fe_values_face[displacement].value(i, q_point)
5177 *
* fe_values_face.JxW(q_point));
5181 *
cell->get_dof_indices(local_dof_indices);
5182 *
constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_matrix, cell_rhs,
5183 *
local_dof_indices,
5199 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_nonlinear_residual"></a>
5200 * <h4>PlasticityContactProblem::compute_nonlinear_residual</h4>
5204 * The following function computes the nonlinear residual of the equation
5205 * given the current solution (or any other linearization point). This
5206 * is needed in the linear search algorithm where we need to
try various
5207 * linear combinations of previous and current (trial) solution to
5208 * compute the (real, globalized) solution of the current Newton step.
5212 * That said, in a slight abuse of the name of the function, it actually
5213 * does significantly more. For example, it also computes the vector
5214 * that corresponds to the Newton residual but without eliminating
5215 * constrained degrees of freedom. We need
this vector to compute contact
5216 * forces and, ultimately, to compute the next active
set. Likewise, by
5217 * keeping track of how many quadrature points we encounter on each cell
5218 * that show plastic yielding, we also compute the
5219 * <code>fraction_of_plastic_q_points_per_cell</code> vector that we
5220 * can later output to visualize the plastic zone. In both of these cases,
5221 * the results are not necessary as part of the line search, and so we may
5222 * be wasting a small amount of time computing them. At the same time,
this
5223 * information appears as a natural by-product of what we need to
do here
5224 * anyway, and we want to collect it once at the
end of each Newton
5225 * step, so we may as well
do it here.
5229 * The actual implementation of
this function should be rather obvious:
5232 *
template <
int dim>
5234 *
ElastoPlasticProblem<dim>::
5238 *
if (base_mesh ==
"Timoshenko beam")
5240 *
traction_surface_id = 5;
5242 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5244 *
traction_surface_id = 0;
5246 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5248 *
traction_surface_id = 2;
5263 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5264 *
const unsigned int n_q_points = quadrature_formula.size();
5265 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5267 *
const EquationData::BodyForce<dim> body_force;
5268 *
std::vector<Vector<double> > body_force_values(n_q_points,
5271 *
const EquationData::
5272 *
IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5273 *
std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5278 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5282 *
newton_rhs_residual = 0;
5284 *
fraction_of_plastic_q_points_per_cell = 0;
5287 *
cell = dof_handler.begin_active(),
5288 *
endc = dof_handler.end();
5289 *
unsigned int cell_number = 0;
5290 *
for (; cell != endc; ++cell, ++cell_number)
5291 *
if (cell->is_locally_owned())
5293 *
fe_values.reinit(cell);
5296 *
std::vector<SymmetricTensor<2, dim> > strain_tensors(n_q_points);
5297 *
fe_values[displacement].get_function_symmetric_gradients(linearization_point,
5300 *
body_force.vector_value_list(fe_values.get_quadrature_points(),
5301 *
body_force_values);
5303 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5306 *
const bool q_point_is_plastic
5307 *
= constitutive_law.get_stress_strain_tensor(strain_tensors[q_point],
5308 *
stress_strain_tensor);
5309 *
if (q_point_is_plastic)
5310 *
++fraction_of_plastic_q_points_per_cell(cell_number);
5313 *
for (
unsigned int i = 0; i < dim; ++i)
5315 *
rhs_values_body_force[i] = body_force_values[q_point][i];
5318 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5320 *
cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5321 *
* rhs_values_body_force
5323 *
strain_tensors[q_point]
5324 *
* stress_strain_tensor
5325 *
* fe_values[displacement].symmetric_gradient(i, q_point)
5327 *
* fe_values.JxW(q_point);
5331 *
cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5333 *
* fe_values.JxW(q_point));
5337 *
for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5338 *
if (cell->face(face)->at_boundary()
5339 *
&& cell->face(face)->boundary_id() == traction_surface_id)
5341 *
fe_values_face.reinit(cell, face);
5343 *
boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5344 *
boundary_force_values);
5346 *
for (
unsigned int q_point = 0; q_point < n_face_q_points;
5350 *
for (
unsigned int i = 0; i < dim; ++i)
5352 *
rhs_values[i] = boundary_force_values[q_point][i];
5354 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5355 *
cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values
5356 *
* fe_values_face.JxW(q_point));
5360 *
cell->get_dof_indices(local_dof_indices);
5361 *
constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_rhs,
5362 *
local_dof_indices,
5363 *
newton_rhs_residual);
5367 *
fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
5379 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton_system"></a>
5380 * <h4>PlasticityContactProblem::solve_newton_system</h4>
5384 * The last piece before we can discuss the actual Newton iteration
5385 * on a single mesh is the solver
for the linear systems. There are
5386 * a couple of complications that slightly obscure the code, but
5387 * mostly it is just setup then solve. Among the complications are:
5391 * - For the hanging nodes we have to
apply
5393 * This is necessary
if a hanging node with solution
value @f$x_0@f$
5394 * has
one neighbor with
value @f$x_1@f$ which is in contact with the
5395 * obstacle and
one neighbor @f$x_2@f$ which is not in contact. Because
5396 * the update
for the former will be prescribed, the hanging node constraint
5397 * will have an inhomogeneity and will look like @f$x_0 = x_1/2 + \text{gap}/2@f$.
5398 * So the corresponding entries in the
5399 * ride-hang-side are non-
zero with a
5402 * - Like in @ref step_40
"step-40", we need to shuffle between vectors that
do and
do
5403 *
do not have ghost elements when solving or
using the solution.
5407 * The rest of the function is similar to @ref step_40
"step-40" and
5408 * @ref step_41
"step-41" except that we use a BiCGStab solver
5409 * instead of CG. This is due to the fact that
for very small hardening
5410 * parameters @f$\gamma@f$, the linear system becomes almost semidefinite though
5411 * still
symmetric. BiCGStab appears to have an easier time with such linear
5415 *
template <
int dim>
5417 *
ElastoPlasticProblem<dim>::solve_newton_system ()
5422 *
distributed_solution = incremental_displacement;
5424 *
constraints_hanging_nodes.set_zero(distributed_solution);
5425 *
constraints_hanging_nodes.set_zero(newton_rhs);
5429 * ------- Solver Bicgstab --- Preconditioner AMG -------------------
5436 * std::vector<std::vector<bool> > constant_modes;
5446 * additional_data.
w_cycle =
false;
5453 * preconditioner.initialize(newton_matrix, additional_data);
5468 *
const double relative_accuracy = 1
e-2;
5469 *
const double solver_tolerance = relative_accuracy
5470 * * newton_matrix.residual(tmp, distributed_solution,
5476 * solver_tolerance);
5478 * solver.solve(newton_matrix, distributed_solution,
5479 * newton_rhs, preconditioner);
5483 * pcout <<
" Error: " << solver_control.initial_value()
5484 * <<
" -> " << solver_control.last_value() <<
" in "
5485 * << solver_control.last_step() <<
" Bicgstab iterations."
5491 * ------- Solver CG --- Preconditioner SSOR -------------------
5499 *
preconditioner.initialize(newton_matrix, additional_data);
5509 *
const double relative_accuracy = 1
e-8;
5512 *
const double relative_accuracy = 1
e-2;
5513 *
const double solver_tolerance = relative_accuracy
5514 *
* newton_matrix.residual(tmp, distributed_solution,
5520 * solver_tolerance);
5524 *
solver_tolerance);
5526 *
solver.solve(newton_matrix, distributed_solution,
5527 *
newton_rhs, preconditioner);
5529 *
pcout <<
" Error: " << solver_control.initial_value()
5530 *
<<
" -> " << solver_control.last_value() <<
" in "
5531 *
<< solver_control.last_step() <<
" CG iterations."
5536 * ........................................................
5542 *
constraints_dirichlet_and_hanging_nodes.distribute(distributed_solution);
5544 *
incremental_displacement = distributed_solution;
5551 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton"></a>
5552 * <h4>PlasticityContactProblem::solve_newton</h4>
5556 * This is,
finally, the function that implements the damped Newton method
5557 * on the current mesh. There are two nested loops: the outer
loop for the Newton
5558 * iteration and the inner
loop for the line search which
5559 * will be used only
if necessary. To obtain a good and reasonable
5560 * starting
value we solve an elastic problem in very
first Newton step on each
5561 * mesh (or only on the
first mesh
if we transfer solutions between meshes). We
5562 *
do so by setting the yield stress to an unreasonably large
value in these
5563 * iterations and then setting it back to the correct
value in subsequent
5568 * Other than
this, the top part of
this function should be reasonably
5572 *
template <
int dim>
5574 *
ElastoPlasticProblem<dim>::solve_newton ()
5583 *
double residual_norm;
5584 *
double previous_residual_norm = -std::numeric_limits<double>::max();
5587 *
previous_disp_norm = 0;
5589 *
const double correct_sigma = sigma_0;
5591 *
const unsigned int max_newton_iter = 100;
5593 *
for (
unsigned int newton_step = 1; newton_step <= max_newton_iter; ++newton_step)
5595 *
if (newton_step == 1
5597 *
((transfer_solution && timestep_no == 1)
5599 *
!transfer_solution))
5600 *
constitutive_law.set_sigma_0(1e+10);
5602 *
constitutive_law.set_sigma_0(correct_sigma);
5604 *
pcout <<
" " << std::endl;
5605 *
pcout <<
" Newton iteration " << newton_step << std::endl;
5607 *
pcout <<
" Assembling system... " << std::endl;
5608 *
newton_matrix = 0;
5610 *
newton_rhs_residual = 0;
5612 *
tmp_solution = solution;
5613 *
tmp_solution += incremental_displacement;
5614 *
assemble_newton_system(tmp_solution,
5615 *
incremental_displacement);
5617 *
pcout <<
" Solving system... " << std::endl;
5618 *
solve_newton_system();
5622 * It gets a bit more hairy after we have computed the
5623 * trial solution @f$\tilde{\mathbf u}@f$ of the current Newton step.
5624 * We handle a highly nonlinear problem so we have to damp
5625 * Newton
's method using a line search. To understand how we do this,
5626 * recall that in our formulation, we compute a trial solution
5627 * in each Newton step and not the update between old and new solution.
5628 * Since the solution set is a convex set, we will use a line
5629 * search that tries linear combinations of the
5630 * previous and the trial solution to guarantee that the
5631 * damped solution is in our solution set again.
5632 * At most we apply 5 damping steps.
5636 * There are exceptions to when we use a line search. First,
5637 * if this is the first Newton step on any mesh, then we don't have
5638 * any
point to compare the residual to, so we
always accept a full
5639 * step. Likewise,
if this is the
second Newton step on the
first mesh (or
5640 * the
second on any mesh
if we don
't transfer solutions from
5641 * mesh to mesh), then we have computed the first of these steps using
5642 * just an elastic model (see how we set the yield stress sigma to
5643 * an unreasonably large value above). In this case, the first Newton
5644 * solution was a purely elastic one, the second one a plastic one,
5645 * and any linear combination would not necessarily be expected to
5646 * lie in the feasible set -- so we just accept the solution we just
5651 * In either of these two cases, we bypass the line search and just
5652 * update residual and other vectors as necessary.
5655 * if ((newton_step==1)
5657 * (transfer_solution && newton_step == 2 && current_refinement_cycle == 0)
5659 * (!transfer_solution && newton_step == 2))
5661 * tmp_solution = solution;
5662 * tmp_solution += incremental_displacement;
5663 * compute_nonlinear_residual(tmp_solution);
5664 * old_solution = incremental_displacement;
5666 * residual = newton_rhs_residual;
5668 * residual.compress(VectorOperation::insert);
5670 * residual_norm = residual.l2_norm();
5672 * pcout << " Accepting Newton solution with residual: "
5673 * << residual_norm << std::endl;
5677 * for (unsigned int i = 0; i < 5; ++i)
5679 * distributed_solution = incremental_displacement;
5681 * const double alpha = std::pow(0.5, static_cast<double>(i));
5682 * tmp_vector = old_solution;
5683 * tmp_vector.sadd(1 - alpha, alpha, distributed_solution);
5685 * TimerOutput::Scope t(computing_timer, "Residual and lambda");
5687 * locally_relevant_tmp_vector = tmp_vector;
5688 * tmp_solution = solution;
5689 * tmp_solution += locally_relevant_tmp_vector;
5690 * compute_nonlinear_residual(tmp_solution);
5691 * residual = newton_rhs_residual;
5693 * residual.compress(VectorOperation::insert);
5695 * residual_norm = residual.l2_norm();
5697 * pcout << " Residual of the system: "
5698 * << residual_norm << std::endl
5699 * << " with a damping parameter alpha = " << alpha
5702 * if (residual_norm < previous_residual_norm)
5706 * incremental_displacement = tmp_vector;
5707 * old_solution = incremental_displacement;
5710 * disp_norm = incremental_displacement.l2_norm();
5715 * The final step is to check for convergence. If the residual is
5716 * less than a threshold of @f$10^{-10}@f$, then we terminate
5717 * the iteration on the current mesh:
5718 * if (residual_norm < 1e-10)
5721 * if (residual_norm < 1e-7)
5724 * pcout << " difference of two consecutive incremental displacement l2 norm : "
5725 * << std::abs(disp_norm - previous_disp_norm) << std::endl;
5726 * if ( std::abs(disp_norm - previous_disp_norm) < 1e-10 &&
5727 * (residual_norm < 1e-5 || std::abs(residual_norm - previous_residual_norm)<1e-9) )
5729 * pcout << " Convergence by difference of two consecutive solution! " << std::endl;
5734 * previous_residual_norm = residual_norm;
5735 * previous_disp_norm = disp_norm;
5742 * <a name="elastoplastic.cc-PlasticityContactProblemcompute_error"></a>
5743 * <h4>PlasticityContactProblem::compute_error</h4>
5749 * template <int dim>
5751 * ElastoPlasticProblem<dim>::compute_error ()
5753 * TrilinosWrappers::MPI::Vector tmp_solution(locally_owned_dofs, mpi_communicator);
5754 * tmp_solution = solution;
5755 * tmp_solution += incremental_displacement;
5757 * estimated_error_per_cell.reinit (triangulation.n_active_cells());
5758 * if (error_estimation_strategy == ErrorEstimationStrategy::kelly_error)
5760 * using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
5762 * KellyErrorEstimator<dim>::estimate(dof_handler,
5763 * QGauss<dim - 1>(fe.degree + 2),
5764 * std::map<types::boundary_id, const Function<dim> *>(),
5766 * estimated_error_per_cell);
5769 * else if (error_estimation_strategy == ErrorEstimationStrategy::residual_error)
5771 * compute_error_residual(tmp_solution);
5774 * else if (error_estimation_strategy == ErrorEstimationStrategy::weighted_residual_error)
5778 * make a non-parallel copy of tmp_solution
5781 * Vector<double> copy_solution(tmp_solution);
5785 * the dual function definition (it should be defined previously, e.g. input file)
5788 * if (base_mesh == "Timoshenko beam")
5790 * double length = .48,
5793 * const Point<dim> evaluation_point(length, -depth/2);
5795 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5797 * DualSolver<dim> dual_solver(triangulation, fe,
5799 * constitutive_law, dual_functional,
5800 * timestep_no, output_dir, base_mesh,
5801 * present_time, end_time);
5803 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5806 * else if (base_mesh == "Thick_tube_internal_pressure")
5808 * const unsigned int face_id = 0;
5809 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5810 * for (unsigned int i=0; i!=dim; ++i)
5812 * comp_stress[i].resize(dim);
5813 * for (unsigned int j=0; j!=dim; ++j)
5815 * comp_stress[i][j] = 1;
5819 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5821 * DualSolver<dim> dual_solver(triangulation, fe,
5823 * constitutive_law, dual_functional,
5824 * timestep_no, output_dir, base_mesh,
5825 * present_time, end_time);
5827 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5830 * else if (base_mesh == "Perforated_strip_tension")
5834 * .........................................
5835 * Mean stress_yy over the bottom boundary
5838 * const unsigned int face_id = 1;
5839 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5840 * for (unsigned int i=0; i!=dim; ++i)
5842 * comp_stress[i].resize(dim);
5843 * for (unsigned int j=0; j!=dim; ++j)
5845 * comp_stress[i][j] = 0;
5848 * comp_stress[1][1] = 1;
5850 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5854 * .........................................
5860 * DualSolver<dim> dual_solver(triangulation, fe,
5862 * constitutive_law, dual_functional,
5863 * timestep_no, output_dir, base_mesh,
5864 * present_time, end_time);
5866 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5869 * else if (base_mesh == "Cantiliver_beam_3d")
5873 * Quantity of interest:
5874 * -----------------------------------------------------------
5875 * displacement at Point A (x=0, y=height/2, z=length)
5879 * const double length = .7,
5882 * const Point<dim> evaluation_point(0, height/2, length);
5884 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5889 * -----------------------------------------------------------
5890 * Mean stress at the specified domain is of interest.
5891 * The interest domains are located on the bottom and top of the flanges
5892 * close to the clamped face, z = 0
5893 * top domain: height/2 - thickness_flange <= y <= height/2
5894 * 0 <= z <= 2 * thickness_flange
5895 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
5896 * 0 <= z <= 2 * thickness_flange
5902 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5903 * for (unsigned int i=0; i!=dim; ++i)
5905 * comp_stress[i].resize(dim);
5906 * for (unsigned int j=0; j!=dim; ++j)
5908 * comp_stress[i][j] = 1;
5911 * DualFunctional::MeanStressDomain<dim> dual_functional(base_mesh, comp_stress);
5915 * -----------------------------------------------------------
5921 * DualSolver<dim> dual_solver(triangulation, fe,
5923 * constitutive_law, dual_functional,
5924 * timestep_no, output_dir, base_mesh,
5925 * present_time, end_time);
5927 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5932 * AssertThrow(false, ExcNotImplemented());
5939 * AssertThrow(false, ExcNotImplemented());
5943 * relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm();
5945 * pcout << "Estimated relative error = " << relative_error << std::endl;
5949 * template <int dim>
5951 * ElastoPlasticProblem<dim>::compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution)
5953 * FEValues<dim> fe_values(fe, quadrature_formula,
5955 * update_gradients |
5957 * update_quadrature_points |
5958 * update_JxW_values);
5960 * const unsigned int n_q_points = quadrature_formula.size();
5961 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5962 * SymmetricTensor<4, dim> stress_strain_tensor_linearized;
5963 * SymmetricTensor<4, dim> stress_strain_tensor;
5964 * Tensor<5, dim> stress_strain_tensor_grad;
5965 * std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
5966 * for (unsigned int i=0; i!=n_q_points; ++i)
5968 * cell_hessians[i].resize (dim);
5970 * const EquationData::BodyForce<dim> body_force;
5972 * std::vector<Vector<double> > body_force_values (n_q_points, Vector<double>(dim));
5973 * const FEValuesExtractors::Vector displacement(0);
5976 * FEFaceValues<dim> fe_face_values_cell(fe, face_quadrature_formula,
5978 * update_quadrature_points|
5979 * update_gradients |
5980 * update_JxW_values |
5981 * update_normal_vectors),
5982 * fe_face_values_neighbor (fe, face_quadrature_formula,
5984 * update_gradients |
5985 * update_JxW_values |
5986 * update_normal_vectors);
5987 * FESubfaceValues<dim> fe_subface_values_cell (fe, face_quadrature_formula,
5988 * update_gradients);
5990 * const unsigned int n_face_q_points = face_quadrature_formula.size();
5991 * std::vector<Vector<double> > jump_residual (n_face_q_points, Vector<double>(dim));
5992 * std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
5993 * for (unsigned int i=0; i!=n_face_q_points; ++i)
5995 * cell_grads[i].resize (dim);
5997 * std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
5998 * for (unsigned int i=0; i!=n_face_q_points; ++i)
6000 * neighbor_grads[i].resize (dim);
6002 * SymmetricTensor<2, dim> q_cell_strain_tensor;
6003 * SymmetricTensor<2, dim> q_neighbor_strain_tensor;
6004 * SymmetricTensor<4, dim> cell_stress_strain_tensor;
6005 * SymmetricTensor<4, dim> neighbor_stress_strain_tensor;
6008 * typename std::map<typename DoFHandler<dim>::face_iterator, Vector<double> >
6010 * typename DoFHandler<dim>::active_cell_iterator
6011 * cell = dof_handler.begin_active(),
6012 * endc = dof_handler.end();
6013 * for (; cell!=endc; ++cell)
6014 * if (cell->is_locally_owned())
6016 * for (unsigned int face_no=0;
6017 * face_no<GeometryInfo<dim>::faces_per_cell;
6020 * face_integrals[cell->face(face_no)].reinit (dim);
6021 * face_integrals[cell->face(face_no)] = -1e20;
6025 * std::vector<Vector<float> > error_indicators_vector;
6026 * error_indicators_vector.resize( triangulation.n_active_cells(),
6027 * Vector<float>(dim) );
6031 * ----------------- estimate_some -------------------------
6034 * cell = dof_handler.begin_active();
6035 * unsigned int present_cell = 0;
6036 * for (; cell!=endc; ++cell, ++present_cell)
6037 * if (cell->is_locally_owned())
6041 * --------------- integrate_over_cell -------------------
6044 * fe_values.reinit(cell);
6045 * body_force.vector_value_list(fe_values.get_quadrature_points(),
6046 * body_force_values);
6047 * fe_values[displacement].get_function_symmetric_gradients(tmp_solution,
6049 * fe_values.get_function_hessians(tmp_solution, cell_hessians);
6051 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
6053 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
6054 * stress_strain_tensor_linearized,
6055 * stress_strain_tensor);
6056 * constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
6057 * cell_hessians[q_point],
6058 * stress_strain_tensor_grad);
6060 * for (unsigned int i=0; i!=dim; ++i)
6062 * error_indicators_vector[present_cell](i) +=
6063 * body_force_values[q_point](i)*fe_values.JxW(q_point);
6064 * for (unsigned int j=0; j!=dim; ++j)
6066 * for (unsigned int k=0; k!=dim; ++k)
6068 * for (unsigned int l=0; l!=dim; ++l)
6070 * error_indicators_vector[present_cell](i) +=
6071 * ( stress_strain_tensor[i][j][k][l]*
6072 * 0.5*(cell_hessians[q_point][k][l][j]
6074 * cell_hessians[q_point][l][k][j])
6075 * + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l]
6077 * fe_values.JxW(q_point);
6087 * -------------------------------------------------------
6088 * compute face_integrals
6091 * for (unsigned int face_no=0;
6092 * face_no<GeometryInfo<dim>::faces_per_cell;
6095 * if (cell->face(face_no)->at_boundary())
6097 * for (unsigned int id=0; id!=dim; ++id)
6099 * face_integrals[cell->face(face_no)](id) = 0;
6104 * if ((cell->neighbor(face_no)->has_children() == false) &&
6105 * (cell->neighbor(face_no)->level() == cell->level()) &&
6106 * (cell->neighbor(face_no)->index() < cell->index()))
6109 * if (cell->at_boundary(face_no) == false)
6110 * if (cell->neighbor(face_no)->level() < cell->level())
6114 * if (cell->face(face_no)->has_children() == false)
6118 * ------------- integrate_over_regular_face -----------
6121 * fe_face_values_cell.reinit(cell, face_no);
6122 * fe_face_values_cell.get_function_gradients (tmp_solution,
6125 * Assert (cell->neighbor(face_no).state() == IteratorState::valid,
6126 * ExcInternalError());
6127 * const unsigned int
6128 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6129 * const typename DoFHandler<dim>::active_cell_iterator
6130 * neighbor = cell->neighbor(face_no);
6132 * fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
6133 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6136 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6138 * q_cell_strain_tensor = 0.;
6139 * q_neighbor_strain_tensor = 0.;
6140 * for (unsigned int i=0; i!=dim; ++i)
6142 * for (unsigned int j=0; j!=dim; ++j)
6144 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6145 * cell_grads[q_point][j][i] );
6146 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6147 * neighbor_grads[q_point][j][i] );
6151 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6152 * cell_stress_strain_tensor);
6153 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6154 * neighbor_stress_strain_tensor);
6156 * jump_residual[q_point] = 0.;
6157 * for (unsigned int i=0; i!=dim; ++i)
6159 * for (unsigned int j=0; j!=dim; ++j)
6161 * for (unsigned int k=0; k!=dim; ++k)
6163 * for (unsigned int l=0; l!=dim; ++l)
6165 * jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
6166 * q_cell_strain_tensor[k][l]
6168 * neighbor_stress_strain_tensor[i][j][k][l]*
6169 * q_neighbor_strain_tensor[k][l] )*
6170 * fe_face_values_cell.normal_vector(q_point)[j];
6178 * Vector<double> face_integral_vector(dim);
6179 * face_integral_vector = 0;
6180 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6182 * for (unsigned int i=0; i!=dim; ++i)
6184 * face_integral_vector(i) += jump_residual[q_point](i) *
6185 * fe_face_values_cell.JxW(q_point);
6189 * Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
6190 * ExcInternalError());
6192 * for (unsigned int i=0; i!=dim; ++i)
6194 * Assert (face_integrals[cell->face(face_no)](i) == -1e20,
6195 * ExcInternalError());
6196 * face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
6202 * -----------------------------------------------------
6210 * ------------- integrate_over_irregular_face ---------
6213 * const typename DoFHandler<dim>::face_iterator
6214 * face = cell->face(face_no);
6215 * const typename DoFHandler<dim>::cell_iterator
6216 * neighbor = cell->neighbor(face_no);
6217 * Assert (neighbor.state() == IteratorState::valid,
6218 * ExcInternalError());
6219 * Assert (neighbor->has_children(),
6220 * ExcInternalError());
6222 * const unsigned int
6223 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6225 * for (unsigned int subface_no=0;
6226 * subface_no<face->n_children(); ++subface_no)
6228 * const typename DoFHandler<dim>::active_cell_iterator
6229 * neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
6230 * Assert (neighbor_child->face(neighbor_neighbor) ==
6231 * cell->face(face_no)->child(subface_no),
6232 * ExcInternalError());
6234 * fe_subface_values_cell.reinit (cell, face_no, subface_no);
6235 * fe_subface_values_cell.get_function_gradients (tmp_solution,
6237 * fe_face_values_neighbor.reinit (neighbor_child,
6238 * neighbor_neighbor);
6239 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6242 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6244 * q_cell_strain_tensor = 0.;
6245 * q_neighbor_strain_tensor = 0.;
6246 * for (unsigned int i=0; i!=dim; ++i)
6248 * for (unsigned int j=0; j!=dim; ++j)
6250 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6251 * cell_grads[q_point][j][i] );
6252 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6253 * neighbor_grads[q_point][j][i] );
6257 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6258 * cell_stress_strain_tensor);
6259 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6260 * neighbor_stress_strain_tensor);
6262 * jump_residual[q_point] = 0.;
6263 * for (unsigned int i=0; i!=dim; ++i)
6265 * for (unsigned int j=0; j!=dim; ++j)
6267 * for (unsigned int k=0; k!=dim; ++k)
6269 * for (unsigned int l=0; l!=dim; ++l)
6271 * jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
6272 * q_cell_strain_tensor[k][l]
6274 * neighbor_stress_strain_tensor[i][j][k][l]*
6275 * q_neighbor_strain_tensor[k][l] )*
6276 * fe_face_values_neighbor.normal_vector(q_point)[j];
6284 * Vector<double> face_integral_vector(dim);
6285 * face_integral_vector = 0;
6286 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6288 * for (unsigned int i=0; i!=dim; ++i)
6290 * face_integral_vector(i) += jump_residual[q_point](i) *
6291 * fe_face_values_neighbor.JxW(q_point);
6295 * for (unsigned int i=0; i!=dim; ++i)
6297 * face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
6302 * Vector<double> sum (dim);
6304 * for (unsigned int subface_no=0;
6305 * subface_no<face->n_children(); ++subface_no)
6307 * Assert (face_integrals.find(face->child(subface_no)) !=
6308 * face_integrals.end(),
6309 * ExcInternalError());
6310 * for (unsigned int i=0; i!=dim; ++i)
6312 * Assert (face_integrals[face->child(subface_no)](i) != -1e20,
6313 * ExcInternalError());
6314 * sum(i) += face_integrals[face->child(subface_no)](i);
6317 * for (unsigned int i=0; i!=dim; ++i)
6319 * face_integrals[face](i) = sum(i);
6325 * -----------------------------------------------------
6335 * ----------------------------------------------------------
6342 * cell = dof_handler.begin_active();
6343 * for (; cell!=endc; ++cell, ++present_cell)
6344 * if (cell->is_locally_owned())
6346 * for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
6349 * Assert(face_integrals.find(cell->face(face_no)) !=
6350 * face_integrals.end(),
6351 * ExcInternalError());
6353 * for (unsigned int id=0; id!=dim; ++id)
6355 * error_indicators_vector[present_cell](id)
6356 * -= 0.5*face_integrals[cell->face(face_no)](id);
6361 * estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
6371 * <a name="elastoplastic.cc-PlasticityContactProblemrefine_grid"></a>
6372 * <h4>PlasticityContactProblem::refine_grid</h4>
6376 * If you've made it
this far into the deal.II tutorial, the following
6377 * function refining the mesh should not pose any challenges to you
6378 * any more. It refines the mesh, either globally or
using the Kelly
6379 * error estimator, and
if so asked also transfers the solution from
6380 * the previous to the next mesh. In the latter
case, we also need
6381 * to compute the active
set and other quantities again,
for which we
6382 * need the information computed by <code>compute_nonlinear_residual()</code>.
6385 *
template <int dim>
6387 *
ElastoPlasticProblem<dim>::refine_grid ()
6391 * ---------------------------------------------------------------
6392 * Make a field variable
for history variables to be able to
6393 * transfer the
data to the quadrature points of the
new mesh
6398 *
history_dof_handler.distribute_dofs (history_fe);
6399 *
std::vector< std::vector< Vector<double> > >
6400 *
history_stress_field (dim, std::vector<
Vector<double> >(dim)),
6401 *
local_history_stress_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6402 *
local_history_stress_fe_values (dim, std::vector<
Vector<double> >(dim));
6405 *
std::vector< std::vector< Vector<double> > >
6406 *
history_strain_field (dim, std::vector<
Vector<double> >(dim)),
6407 *
local_history_strain_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6408 *
local_history_strain_fe_values (dim, std::vector<
Vector<double> >(dim));
6410 *
for (
unsigned int i=0; i<dim; ++i)
6411 *
for (
unsigned int j=0; j<dim; ++j)
6413 *
history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6414 *
local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6415 *
local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6417 *
history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6418 *
local_history_strain_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6419 *
local_history_strain_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6422 *
quadrature_formula.size());
6425 *
quadrature_formula, quadrature_formula,
6426 *
qpoint_to_dof_matrix);
6428 *
cell = dof_handler.begin_active(),
6429 *
endc = dof_handler.end(),
6430 *
dg_cell = history_dof_handler.begin_active();
6431 *
for (; cell!=endc; ++cell, ++dg_cell)
6432 *
if (cell->is_locally_owned())
6434 *
PointHistory<dim> *local_quadrature_points_history
6435 *
=
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6436 *
Assert (local_quadrature_points_history >=
6437 *
&quadrature_point_history.front(),
6438 *
ExcInternalError());
6439 *
Assert (local_quadrature_points_history <
6440 *
&quadrature_point_history.back(),
6441 *
ExcInternalError());
6442 *
for (
unsigned int i=0; i<dim; ++i)
6443 *
for (
unsigned int j=0; j<dim; ++j)
6445 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6447 *
local_history_stress_values_at_qpoints[i][j](q)
6448 *
= local_quadrature_points_history[q].old_stress[i][j];
6450 *
local_history_strain_values_at_qpoints[i][j](q)
6451 *
= local_quadrature_points_history[q].old_strain[i][j];
6453 *
qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
6454 *
local_history_stress_values_at_qpoints[i][j]);
6455 *
dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
6456 *
history_stress_field[i][j]);
6458 *
qpoint_to_dof_matrix.vmult (local_history_strain_fe_values[i][j],
6459 *
local_history_strain_values_at_qpoints[i][j]);
6460 *
dg_cell->set_dof_values (local_history_strain_fe_values[i][j],
6461 *
history_strain_field[i][j]);
6468 * ---------------------------------------------------------------
6472 *
if (refinement_strategy == RefinementStrategy::refine_global)
6475 *
cell = triangulation.begin_active();
6476 *
cell != triangulation.end(); ++cell)
6477 *
if (cell->is_locally_owned())
6478 *
cell->set_refine_flag ();
6482 *
const double refine_fraction_cells = .3,
6483 *
coarsen_fraction_cells = .03;
6486 *
const double refine_fraction_cells = .1,
6487 * coarsen_fraction_cells = .3;
6494 * ::refine_and_coarsen_fixed_number(triangulation,
6495 *
estimated_error_per_cell,
6496 *
refine_fraction_cells, coarsen_fraction_cells);
6499 *
triangulation.prepare_coarsening_and_refinement();
6502 *
solution_transfer.prepare_for_coarsening_and_refinement(solution);
6505 *
if (transfer_solution)
6506 *
incremental_displacement_transfer.prepare_for_coarsening_and_refinement(incremental_displacement);
6509 *
history_stress_field_transfer1(history_dof_handler),
6510 *
history_stress_field_transfer2(history_dof_handler);
6511 *
history_stress_field_transfer0.prepare_for_coarsening_and_refinement(history_stress_field[0]);
6514 *
history_stress_field_transfer1.prepare_for_coarsening_and_refinement(history_stress_field[1]);
6518 *
history_stress_field_transfer2.prepare_for_coarsening_and_refinement(history_stress_field[2]);
6522 *
history_strain_field_transfer1(history_dof_handler),
6523 *
history_strain_field_transfer2(history_dof_handler);
6524 *
history_strain_field_transfer0.prepare_for_coarsening_and_refinement(history_strain_field[0]);
6527 *
history_strain_field_transfer1.prepare_for_coarsening_and_refinement(history_strain_field[1]);
6531 *
history_strain_field_transfer2.prepare_for_coarsening_and_refinement(history_strain_field[2]);
6534 *
triangulation.execute_coarsening_and_refinement();
6535 *
pcout <<
" Number of active cells: "
6536 *
<< triangulation.n_active_cells()
6540 *
setup_quadrature_point_history ();
6546 * distributed_solution = solution;
6549 *
solution_transfer.interpolate(distributed_solution);
6550 *
solution = distributed_solution;
6552 *
if (transfer_solution)
6557 * distributed_incremental_displacement = incremental_displacement;
6560 *
incremental_displacement_transfer.interpolate(distributed_incremental_displacement);
6561 *
incremental_displacement = distributed_incremental_displacement;
6564 * compute_nonlinear_residual(incremental_displacement);
6571 * ---------------------------------------------------
6574 *
history_dof_handler.distribute_dofs (history_fe);
6581 *
for (
unsigned int i=0; i<dim; ++i)
6582 *
for (
unsigned int j=0; j<dim; ++j)
6584 *
history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6587 *
history_stress_field_transfer0.interpolate(history_stress_field[0]);
6590 *
history_stress_field_transfer1.interpolate(history_stress_field[1]);
6594 *
history_stress_field_transfer2.interpolate(history_stress_field[2]);
6602 *
for (
unsigned int i=0; i<dim; ++i)
6603 *
for (
unsigned int j=0; j<dim; ++j)
6605 *
history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6608 *
history_strain_field_transfer0.interpolate(history_strain_field[0]);
6611 *
history_strain_field_transfer1.interpolate(history_strain_field[1]);
6615 *
history_strain_field_transfer2.interpolate(history_strain_field[2]);
6620 * ---------------------------------------------------------------
6621 * Transfer the history
data to the quadrature points of the
new mesh
6622 * In a
final step, we have to get the
data back from the now
6623 * interpolated global field to the quadrature points on the
6624 *
new mesh. The following code will
do that:
6631 *
history_fe.dofs_per_cell);
6634 *
quadrature_formula,
6635 *
dof_to_qpoint_matrix);
6636 *
cell = dof_handler.begin_active();
6637 *
endc = dof_handler.end();
6638 *
dg_cell = history_dof_handler.begin_active();
6639 *
for (; cell != endc; ++cell, ++dg_cell)
6640 *
if (cell->is_locally_owned())
6642 *
PointHistory<dim> *local_quadrature_points_history
6643 *
=
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6644 *
Assert (local_quadrature_points_history >=
6645 *
&quadrature_point_history.front(),
6646 *
ExcInternalError());
6647 *
Assert (local_quadrature_points_history <
6648 *
&quadrature_point_history.back(),
6649 *
ExcInternalError());
6650 *
for (
unsigned int i=0; i<dim; ++i)
6651 *
for (
unsigned int j=0; j<dim; ++j)
6653 *
dg_cell->get_dof_values (history_stress_field[i][j],
6654 *
local_history_stress_fe_values[i][j]);
6655 *
dof_to_qpoint_matrix.vmult (local_history_stress_values_at_qpoints[i][j],
6656 *
local_history_stress_fe_values[i][j]);
6658 *
dg_cell->get_dof_values (history_strain_field[i][j],
6659 *
local_history_strain_fe_values[i][j]);
6660 *
dof_to_qpoint_matrix.vmult (local_history_strain_values_at_qpoints[i][j],
6661 *
local_history_strain_fe_values[i][j]);
6662 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6664 *
local_quadrature_points_history[q].old_stress[i][j]
6665 *
= local_history_stress_values_at_qpoints[i][j](q);
6667 *
local_quadrature_points_history[q].old_strain[i][j]
6668 *
= local_history_strain_values_at_qpoints[i][j](q);
6679 * <a name=
"elastoplastic.cc-ElastoPlasticProblemsetup_quadrature_point_history"></a>
6680 * <h4>ElastoPlasticProblem::setup_quadrature_point_history</h4>
6685 * of the history variables, such as the existing stresses in the material,
6686 * that we store in each quadrature
point. As mentioned above, we use the
6687 * <code>user_pointer</code>
for this that is available in each cell.
6691 * To put
this into larger perspective, we note that
if we had previously
6692 * available stresses in our model (which we assume
do not exist
for the
6693 * purpose of
this program), then we would need to
interpolate the field of
6694 * preexisting stresses to the quadrature points. Likewise,
if we were to
6695 * simulate elasto-plastic materials with hardening/softening, then we would
6696 * have to store additional history variables like the present yield stress
6697 * of the accumulated plastic strains in each quadrature
6698 * points. Pre-existing hardening or weakening would then be implemented by
6699 * interpolating these variables in the present function as well.
6702 *
template <
int dim>
6703 *
void ElastoPlasticProblem<dim>::setup_quadrature_point_history ()
6707 * What we need to
do here is to
first count how many quadrature points
6708 * are within the responsibility of
this processor. This, of course,
6709 * equals the number of cells that belong to
this processor times the
6710 * number of quadrature points our quadrature formula has on each cell.
6714 * For good measure, we also
set all user pointers of all cells, whether
6715 * ours of not, to the null
pointer. This way,
if we ever access the user
6716 *
pointer of a cell which we should not have accessed, a segmentation
6717 * fault will let us know that
this should not have happened:
6720 *
unsigned int our_cells = 0;
6722 *
cell = triangulation.begin_active();
6723 *
cell != triangulation.end(); ++cell)
6724 *
if (cell->is_locally_owned())
6727 *
triangulation.clear_user_data();
6731 * Next, allocate as many quadrature objects as we need. Since the
6732 * <code>resize</code> function does not actually shrink the amount of
6733 * allocated memory
if the requested
new size is smaller than the old
6734 *
size, we resort to a trick to
first free all memory, and then
6735 * reallocate it: we declare an empty vector as a temporary variable and
6736 * then
swap the contents of the old vector and
this temporary
6737 * variable. This makes sure that the
6738 * <code>quadrature_point_history</code> is now really empty, and we can
6739 * let the temporary variable that now holds the previous contents of the
6740 * vector go out of scope and be destroyed. In the next step. we can then
6741 * re-allocate as many elements as we need, with the vector
6742 *
default-initializing the <code>PointHistory</code> objects, which
6743 * includes setting the stress variables to
zero.
6747 *
std::vector<PointHistory<dim> > tmp;
6748 *
tmp.swap (quadrature_point_history);
6750 *
quadrature_point_history.resize (our_cells *
6751 *
quadrature_formula.size());
6755 * Finally
loop over all cells again and
set the user pointers from the
6756 * cells that belong to the present processor to
point to the
first
6757 * quadrature
point objects corresponding to
this cell in the vector of
6761 *
unsigned int history_index = 0;
6763 *
cell = triangulation.begin_active();
6764 *
cell != triangulation.end(); ++cell)
6765 *
if (cell->is_locally_owned())
6767 *
cell->set_user_pointer (&quadrature_point_history[history_index]);
6768 *
history_index += quadrature_formula.size();
6773 * At the
end,
for good measure make sure that our count of elements was
6774 * correct and that we have both used up all objects we allocated
6775 * previously, and not
point to any objects beyond the
end of the
6776 * vector. Such defensive programming strategies are
always good checks to
6777 * avoid accidental errors and to guard against future changes to
this
6778 * function that forget to update all uses of a variable at the same
6779 * time. Recall that constructs
using the <code>
Assert</code> macro are
6780 * optimized away in optimized mode, so
do not affect the
run time of
6784 *
Assert (history_index == quadrature_point_history.size(),
6785 *
ExcInternalError());
6791 * <a name=
"elastoplastic.cc-ElastoPlasticProblemupdate_quadrature_point_history"></a>
6792 * <h4>ElastoPlasticProblem::update_quadrature_point_history</h4>
6796 * At the
end of each time step, we should have computed an incremental
6797 * displacement update so that the material in its
new configuration
6798 * accommodates
for the difference between the external body and boundary
6799 * forces applied during
this time step minus the forces exerted through
6800 * preexisting
internal stresses. In order to have the preexisting
6801 * stresses available at the next time step, we therefore have to update the
6802 * preexisting stresses with the stresses due to the incremental
6803 * displacement computed during the present time step. Ideally, the
6804 * resulting
sum of
internal stresses would exactly counter all external
6805 * forces. Indeed, a simple experiment can make sure that
this is so:
if we
6806 * choose boundary conditions and body forces to be time
independent, then
6807 * the forcing terms (the sum of external forces and
internal stresses)
6808 * should be exactly
zero. If you make
this experiment, you will realize
6809 * from the output of the
norm of the right hand side in each time step that
6810 *
this is almost the
case: it is not exactly
zero, since in the
first time
6811 * step the incremental displacement and stress updates were computed
6812 * relative to the undeformed mesh, which was then deformed. In the
second
6813 * time step, we again compute displacement and stress updates, but
this
6814 * time in the deformed mesh -- there, the resulting updates are very small
6815 * but not quite
zero. This can be iterated, and in each such iteration the
6816 * residual, i.e. the
norm of the right hand side vector, is reduced;
if one
6817 * makes
this little experiment,
one realizes that the
norm of
this residual
6818 * decays exponentially with the number of iterations, and after an
initial
6819 * very rapid decline is reduced by roughly a factor of about 3.5 in each
6820 * iteration (
for one testcase I looked at, other testcases, and other
6821 *
numbers of unknowns change the factor, but not the exponential decay).
6825 * In a sense,
this can then be considered as a quasi-timestepping scheme to
6826 * resolve the nonlinear problem of solving large-deformation elasticity on
6827 * a mesh that is moved along in a Lagrangian manner.
6831 * Another complication is that the existing (old) stresses are defined on
6832 * the old mesh, which we will move around after updating the stresses. If
6833 *
this mesh update involves rotations of the cell, then we need to also
6834 *
rotate the updated stress, since it was computed relative to the
6835 * coordinate system of the old cell.
6839 * Thus, what we need is the following: on each cell which the present
6840 * processor owns, we need to
extract the old stress from the
data stored
6841 * with each quadrature
point, compute the stress update, add the two
6842 * together, and then
rotate the result together with the incremental
6843 * rotation computed from the incremental displacement at the present
6844 * quadrature
point. We will detail these steps below:
6847 *
template <
int dim>
6848 *
void ElastoPlasticProblem<dim>::
6849 *
update_quadrature_point_history ()
6853 * First,
set up an <code>
FEValues</code>
object by which we will evaluate
6854 * the displacements and the
gradients thereof at the
6855 * quadrature points, together with a vector that will hold
this
6863 *
const unsigned int n_q_points = quadrature_formula.size();
6865 *
std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
6871 * Then
loop over all cells and
do the job in the cells that belong to our
6879 *
cell = dof_handler.begin_active(),
6880 *
endc = dof_handler.end();
6884 *
for (; cell != endc; ++cell)
6885 *
if (cell->is_locally_owned())
6890 * the present cell, and, as a defensive measure, make sure that
6891 *
this pointer is within the bounds of the global array:
6894 *
PointHistory<dim> *local_quadrature_points_history
6895 *
=
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6896 *
Assert (local_quadrature_points_history >=
6897 *
&quadrature_point_history.front(),
6898 *
ExcInternalError());
6899 *
Assert (local_quadrature_points_history <
6900 *
&quadrature_point_history.back(),
6901 *
ExcInternalError());
6905 * Then initialize the <code>
FEValues</code>
object on the present
6906 * cell, and
extract the strains of the displacement at the
6910 *
fe_values.reinit (cell);
6911 *
fe_values[displacement].get_function_symmetric_gradients(incremental_displacement,
6912 *
incremental_strain_tensor);
6916 * Then
loop over the quadrature points of
this cell:
6919 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6921 *
local_quadrature_points_history[q].old_strain +=
6922 *
incremental_strain_tensor[q];
6924 *
constitutive_law.get_stress_strain_tensor(local_quadrature_points_history[q].old_strain,
6925 *
stress_strain_tensor);
6929 * The result of these operations is then written back into
6930 * the original place:
6933 *
local_quadrature_points_history[q].old_stress
6934 *
= stress_strain_tensor * local_quadrature_points_history[q].old_strain;
6936 *
local_quadrature_points_history[q].point
6937 *
= fe_values.get_quadrature_points ()[q];
6946 * <a name=
"elastoplastic.cc-PlasticityContactProblemmove_mesh"></a>
6947 * <h4>PlasticityContactProblem::move_mesh</h4>
6951 * The remaining three
functions before we get to <code>
run()</code>
6952 * have to
do with generating output. The following
one is an attempt
6953 * at showing the deformed body in its deformed configuration. To
this
6954 *
end,
this function takes a displacement vector field and moves every
6955 * vertex of the (local part) of the mesh by the previously computed
6956 * displacement. We will call
this function with the current
6957 * displacement field before we generate graphical output, and we will
6958 * call it again after generating graphical output with the
negative
6959 * displacement field to undo the changes to the mesh so made.
6963 * The function itself is pretty straightforward. All we have to
do
6964 * is keep track which vertices we have already touched, as we
6965 * encounter the same vertices multiple times as we
loop over cells.
6968 *
template <
int dim>
6970 *
ElastoPlasticProblem<dim>::
6973 *
std::vector<bool> vertex_touched(triangulation.n_vertices(),
false);
6976 *
dof_handler.begin_active();
6977 *
cell != dof_handler.end(); ++cell)
6978 *
if (cell->is_locally_owned())
6979 *
for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
6980 *
if (vertex_touched[cell->vertex_index(v)] ==
false)
6982 *
vertex_touched[cell->vertex_index(v)] =
true;
6985 *
for (
unsigned int d = 0;
d < dim; ++
d)
6986 *
vertex_displacement[d] = displacement(cell->vertex_dof_index(v, d));
6988 *
cell->vertex(v) += vertex_displacement;
6997 * <a name=
"elastoplastic.cc-PlasticityContactProblemoutput_results"></a>
6998 * <h4>PlasticityContactProblem::output_results</h4>
7002 * Next is the function we use to actually generate graphical output. The
7003 * function is a bit tedious, but not actually particularly complicated.
7004 * It moves the mesh at the top (and moves it back at the
end), then
7005 * computes the contact forces along the contact surface. We can
do
7006 * so (as shown in the accompanying paper) by taking the untreated
7007 * residual vector and identifying which degrees of freedom
7008 * correspond to those with contact by asking whether they have an
7009 * inhomogeneous constraints associated with them. As
always, we need
7010 * to be mindful that we can only write into completely distributed
7011 * vectors (i.e., vectors without ghost elements) but that when we
7012 * want to generate output, we need vectors that
do indeed have
7013 * ghost entries
for all locally relevant degrees of freedom.
7016 *
template <
int dim>
7018 *
ElastoPlasticProblem<dim>::output_results (
const std::string &filename_base)
7022 *
pcout <<
" Writing graphical output... " << std::flush;
7026 *
const double magnified_factor = 3;
7027 *
magnified_solution *= magnified_factor;
7029 *
move_mesh(magnified_solution);
7036 *
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
7038 *
data_out.add_data_vector(solution,
7039 *
std::vector<std::string> (dim,
"displacement"),
7043 *
std::vector<std::string> solution_names;
7048 *
solution_names.push_back (
"displacement");
7051 *
solution_names.push_back (
"x_displacement");
7052 *
solution_names.push_back (
"y_displacement");
7055 *
solution_names.push_back (
"x_displacement");
7056 *
solution_names.push_back (
"y_displacement");
7057 *
solution_names.push_back (
"z_displacement");
7063 *
data_out.add_data_vector (solution, solution_names);
7068 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
7069 *
subdomain(i) = triangulation.locally_owned_subdomain();
7070 *
data_out.add_data_vector(subdomain,
"subdomain");
7073 *
data_out.add_data_vector(fraction_of_plastic_q_points_per_cell,
7074 *
"fraction_of_plastic_q_points");
7077 *
data_out.build_patches();
7081 * In the remainder of the function, we generate
one VTU file on
7082 * every processor, indexed by the subdomain
id of
this processor.
7083 * On the
first processor, we then also create a <code>.pvtu</code>
7084 * file that indexes <i>all</i> of the VTU files so that the entire
7085 *
set of output files can be read at once. These <code>.pvtu</code>
7086 * are used by Paraview to describe an entire
parallel computation
's
7087 * output files. We then do the same again for the competitor of
7088 * Paraview, the Visit visualization program, by creating a matching
7089 * <code>.visit</code> file.
7092 * const std::string filename =
7093 * (output_dir + filename_base + "-"
7094 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7096 * std::ofstream output_vtu((filename + ".vtu").c_str());
7097 * data_out.write_vtu(output_vtu);
7098 * pcout << output_dir + filename_base << ".pvtu" << std::endl;
7101 * if (this_mpi_process == 0)
7103 * std::vector<std::string> filenames;
7104 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7105 * filenames.push_back(filename_base + "-" +
7106 * Utilities::int_to_string(i, 4) +
7109 * std::ofstream pvtu_master_output((output_dir + filename_base + ".pvtu").c_str());
7110 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7112 * std::ofstream visit_master_output((output_dir + filename_base + ".visit").c_str());
7113 * data_out.write_pvtu_record(visit_master_output, filenames);
7117 * produce eps files for mesh illustration
7120 * std::ofstream output_eps((filename + ".eps").c_str());
7122 * grid_out.write_eps(triangulation, output_eps);
7127 * Extrapolate the stresses from Gauss point to the nodes
7130 * SymmetricTensor<2, dim> stress_at_qpoint;
7132 * FE_DGQ<dim> history_fe (1);
7133 * DoFHandler<dim> history_dof_handler (triangulation);
7134 * history_dof_handler.distribute_dofs (history_fe);
7135 * std::vector< std::vector< Vector<double> > >
7136 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7137 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7138 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7139 * for (unsigned int i=0; i<dim; ++i)
7140 * for (unsigned int j=0; j<dim; ++j)
7142 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7143 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7144 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7147 * Vector<double> VM_stress_field (history_dof_handler.n_dofs()),
7148 * local_VM_stress_values_at_qpoints (quadrature_formula.size()),
7149 * local_VM_stress_fe_values (history_fe.dofs_per_cell);
7151 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7152 * quadrature_formula.size());
7153 * FETools::compute_projection_from_quadrature_points_matrix
7155 * quadrature_formula, quadrature_formula,
7156 * qpoint_to_dof_matrix);
7158 * typename DoFHandler<dim>::active_cell_iterator
7159 * cell = dof_handler.begin_active(),
7160 * endc = dof_handler.end(),
7161 * dg_cell = history_dof_handler.begin_active();
7163 * const FEValuesExtractors::Vector displacement(0);
7165 * for (; cell!=endc; ++cell, ++dg_cell)
7166 * if (cell->is_locally_owned())
7168 * PointHistory<dim> *local_quadrature_points_history
7169 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7170 * Assert (local_quadrature_points_history >=
7171 * &quadrature_point_history.front(),
7172 * ExcInternalError());
7173 * Assert (local_quadrature_points_history <
7174 * &quadrature_point_history.back(),
7175 * ExcInternalError());
7179 * Then loop over the quadrature points of this cell:
7182 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7184 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7186 * for (unsigned int i=0; i<dim; ++i)
7187 * for (unsigned int j=i; j<dim; ++j)
7189 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7192 * local_VM_stress_values_at_qpoints(q) = Evaluation::get_von_Mises_stress(stress_at_qpoint);
7197 * for (unsigned int i=0; i<dim; ++i)
7198 * for (unsigned int j=i; j<dim; ++j)
7200 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7201 * local_history_stress_values_at_qpoints[i][j]);
7202 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7203 * history_stress_field[i][j]);
7206 * qpoint_to_dof_matrix.vmult (local_VM_stress_fe_values,
7207 * local_VM_stress_values_at_qpoints);
7208 * dg_cell->set_dof_values (local_VM_stress_fe_values,
7216 * Save stresses on nodes by nodal averaging
7217 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7218 * in order to compute stresses on nodes (by applying nodal averaging)
7219 * Therefore, each vertex has one degree of freedom
7222 * FE_Q<dim> fe_1 (1);
7223 * DoFHandler<dim> dof_handler_1 (triangulation);
7224 * dof_handler_1.distribute_dofs (fe_1);
7226 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7227 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7229 * std::vector< std::vector< Vector<double> > >
7230 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7231 * for (unsigned int i=0; i<dim; ++i)
7232 * for (unsigned int j=0; j<dim; ++j)
7234 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7237 * Vector<double> VM_stress_on_vertices (dof_handler_1.n_dofs()),
7238 * counter_on_vertices (dof_handler_1.n_dofs());
7239 * VM_stress_on_vertices = 0;
7240 * counter_on_vertices = 0;
7242 * cell = dof_handler.begin_active();
7243 * dg_cell = history_dof_handler.begin_active();
7244 * typename DoFHandler<dim>::active_cell_iterator
7245 * cell_1 = dof_handler_1.begin_active();
7246 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7247 * if (cell->is_locally_owned())
7249 * dg_cell->get_dof_values (VM_stress_field,
7250 * local_VM_stress_fe_values);
7252 * for (unsigned int i=0; i<dim; ++i)
7253 * for (unsigned int j=0; j<dim; ++j)
7255 * dg_cell->get_dof_values (history_stress_field[i][j],
7256 * local_history_stress_fe_values[i][j]);
7259 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7261 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7266 * Point<dim> point1, point2;
7267 * point1 = cell_1->vertex(v);
7268 * point2 = dg_cell->vertex(v);
7269 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7276 * counter_on_vertices (dof_1_vertex) += 1;
7278 * VM_stress_on_vertices (dof_1_vertex) += local_VM_stress_fe_values (v);
7280 * for (unsigned int i=0; i<dim; ++i)
7281 * for (unsigned int j=0; j<dim; ++j)
7283 * history_stress_on_vertices[i][j](dof_1_vertex) +=
7284 * local_history_stress_fe_values[i][j](v);
7290 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7292 * VM_stress_on_vertices(id) /= counter_on_vertices(id);
7294 * for (unsigned int i=0; i<dim; ++i)
7295 * for (unsigned int j=0; j<dim; ++j)
7297 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7303 * Save figures of stresses
7306 * if (show_stresses)
7309 * DataOut<dim> data_out;
7310 * data_out.attach_dof_handler (history_dof_handler);
7313 * data_out.add_data_vector (history_stress_field[0][0], "stress_xx");
7314 * data_out.add_data_vector (history_stress_field[1][1], "stress_yy");
7315 * data_out.add_data_vector (history_stress_field[0][1], "stress_xy");
7316 * data_out.add_data_vector (VM_stress_field, "Von_Mises_stress");
7320 * data_out.add_data_vector (history_stress_field[0][2], "stress_xz");
7321 * data_out.add_data_vector (history_stress_field[1][2], "stress_yz");
7322 * data_out.add_data_vector (history_stress_field[2][2], "stress_zz");
7325 * data_out.build_patches ();
7327 * const std::string filename_base_stress = ("stress-" + filename_base);
7329 * const std::string filename =
7330 * (output_dir + filename_base_stress + "-"
7331 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7333 * std::ofstream output_vtu((filename + ".vtu").c_str());
7334 * data_out.write_vtu(output_vtu);
7335 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7337 * if (this_mpi_process == 0)
7339 * std::vector<std::string> filenames;
7340 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7341 * filenames.push_back(filename_base_stress + "-" +
7342 * Utilities::int_to_string(i, 4) +
7345 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7346 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7348 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7349 * data_out.write_pvtu_record(visit_master_output, filenames);
7356 * DataOut<dim> data_out;
7357 * data_out.attach_dof_handler (dof_handler_1);
7360 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_xx_averaged");
7361 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_yy_averaged");
7362 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_xy_averaged");
7363 * data_out.add_data_vector (VM_stress_on_vertices, "Von_Mises_stress_averaged");
7367 * data_out.add_data_vector (history_stress_on_vertices[0][2], "stress_xz_averaged");
7368 * data_out.add_data_vector (history_stress_on_vertices[1][2], "stress_yz_averaged");
7369 * data_out.add_data_vector (history_stress_on_vertices[2][2], "stress_zz_averaged");
7372 * data_out.build_patches ();
7374 * const std::string filename_base_stress = ("averaged-stress-" + filename_base);
7376 * const std::string filename =
7377 * (output_dir + filename_base_stress + "-"
7378 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7380 * std::ofstream output_vtu((filename + ".vtu").c_str());
7381 * data_out.write_vtu(output_vtu);
7382 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7384 * if (this_mpi_process == 0)
7386 * std::vector<std::string> filenames;
7387 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7388 * filenames.push_back(filename_base_stress + "-" +
7389 * Utilities::int_to_string(i, 4) +
7392 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7393 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7395 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7396 * data_out.write_pvtu_record(visit_master_output, filenames);
7403 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7411 * magnified_solution *= -1;
7412 * move_mesh(magnified_solution);
7419 * if (base_mesh == "Timoshenko beam")
7421 * const double length = .48,
7424 * Point<dim> intersted_point(length, -depth/2);
7425 * Point<dim> vertex_displacement;
7426 * bool vertex_found = false;
7428 * for (typename DoFHandler<dim>::active_cell_iterator cell =
7429 * dof_handler.begin_active();
7430 * cell != dof_handler.end(); ++cell)
7431 * if (cell->is_locally_owned() && !vertex_found)
7432 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7433 * if ( std::fabs(cell->vertex(v)[0] - intersted_point[0])<1e-6 &&
7434 * std::fabs(cell->vertex(v)[1] - intersted_point[1])<1e-6)
7436 * vertex_found = true;
7438 * for (unsigned int d = 0; d < dim; ++d)
7439 * vertex_displacement[d] = solution(cell->vertex_dof_index(v, d));
7444 * pcout << " Number of active cells: "
7445 * << triangulation.n_global_active_cells() << std::endl
7446 * << " Number of degrees of freedom: " << dof_handler.n_dofs()
7449 * AssertThrow(vertex_found, ExcInternalError());
7450 * std::cout << "Displacement at the point (" << intersted_point[0]
7451 * << ", " << intersted_point[1] << ") is "
7452 * << "(" << vertex_displacement[0]
7453 * << ", " << vertex_displacement[1] << ").\n";
7455 * Vector<double> vertex_exact_displacement(dim);
7456 * EquationData::IncrementalBoundaryValues<dim> incremental_boundary_values(present_time, end_time);
7457 * incremental_boundary_values.vector_value (intersted_point, vertex_exact_displacement);
7459 * std::cout << "Exact displacement at the point (" << intersted_point[0]
7460 * << ", " << intersted_point[1] << ") is "
7461 * << "(" << vertex_exact_displacement[0]
7462 * << ", " << vertex_exact_displacement[1] << ").\n\n";
7465 * else if (base_mesh == "Thick_tube_internal_pressure")
7467 * const double pressure (0.6*2.4e8),
7468 * inner_radius (.1);
7471 * const double pressure (1.94e8),
7472 * inner_radius (.1);
7480 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7481 * 3d and plane strain
7484 * const double mu (e_modulus / (2 * (1 + nu)));
7486 * const Point<dim> point_A(inner_radius, 0.);
7487 * Vector<double> disp_A(dim);
7491 * make a non-parallel copy of solution
7494 * Vector<double> copy_solution(solution);
7496 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
7498 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
7500 * table_results.add_value("time step", timestep_no);
7501 * table_results.add_value("Cells", triangulation.n_global_active_cells());
7502 * table_results.add_value("DoFs", dof_handler.n_dofs());
7503 * table_results.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7504 * table_results.add_value("4*mu*u_A/(sigma_0*a)", 4*mu*disp_A(0)/(sigma_0*inner_radius));
7508 * Compute stresses in the POLAR coordinates, 1- save it on Gauss points,
7509 * 2- extrapolate them to nodes and taking their avarages (nodal avaraging)
7512 * AssertThrow (dim == 2, ExcNotImplemented());
7516 * we define a rotation matrix to be able to transform the stress
7517 * from the Cartesian coordinate to the polar coordinate
7520 * Tensor<2, dim> rotation_matrix; // [cos sin; -sin cos] , sigma_r = rot * sigma * rot^T
7522 * FEValues<dim> fe_values (fe, quadrature_formula, update_quadrature_points |
7523 * update_values | update_gradients);
7525 * const unsigned int n_q_points = quadrature_formula.size();
7527 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
7528 * SymmetricTensor<4, dim> stress_strain_tensor;
7529 * Tensor<2, dim> stress_at_qpoint;
7531 * FE_DGQ<dim> history_fe (1);
7532 * DoFHandler<dim> history_dof_handler (triangulation);
7533 * history_dof_handler.distribute_dofs (history_fe);
7534 * std::vector< std::vector< Vector<double> > >
7535 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7536 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7537 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7538 * for (unsigned int i=0; i<dim; ++i)
7539 * for (unsigned int j=0; j<dim; ++j)
7541 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7542 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7543 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7546 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7547 * quadrature_formula.size());
7548 * FETools::compute_projection_from_quadrature_points_matrix
7550 * quadrature_formula, quadrature_formula,
7551 * qpoint_to_dof_matrix);
7553 * typename DoFHandler<dim>::active_cell_iterator
7554 * cell = dof_handler.begin_active(),
7555 * endc = dof_handler.end(),
7556 * dg_cell = history_dof_handler.begin_active();
7558 * const FEValuesExtractors::Vector displacement(0);
7560 * for (; cell!=endc; ++cell, ++dg_cell)
7561 * if (cell->is_locally_owned())
7563 * PointHistory<dim> *local_quadrature_points_history
7564 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7565 * Assert (local_quadrature_points_history >=
7566 * &quadrature_point_history.front(),
7567 * ExcInternalError());
7568 * Assert (local_quadrature_points_history <
7569 * &quadrature_point_history.back(),
7570 * ExcInternalError());
7574 * Then loop over the quadrature points of this cell:
7577 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7579 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7583 * transform the stress from the Cartesian coordinate to the polar coordinate
7586 * const Point<dim> point = local_quadrature_points_history[q].point;
7587 * const double theta = std::atan2(point(1),point(0));
7594 * rotation_matrix[0][0] = std::cos(theta);
7595 * rotation_matrix[0][1] = std::sin(theta);
7596 * rotation_matrix[1][0] = -std::sin(theta);
7597 * rotation_matrix[1][1] = std::cos(theta);
7601 * stress in polar coordinate
7604 * stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix);
7606 * for (unsigned int i=0; i<dim; ++i)
7607 * for (unsigned int j=i; j<dim; ++j)
7609 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7615 * for (unsigned int i=0; i<dim; ++i)
7616 * for (unsigned int j=i; j<dim; ++j)
7618 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7619 * local_history_stress_values_at_qpoints[i][j]);
7620 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7621 * history_stress_field[i][j]);
7627 * DataOut<dim> data_out;
7628 * data_out.attach_dof_handler (history_dof_handler);
7631 * data_out.add_data_vector (history_stress_field[0][0], "stress_rr");
7632 * data_out.add_data_vector (history_stress_field[1][1], "stress_tt");
7633 * data_out.add_data_vector (history_stress_field[0][1], "stress_rt");
7635 * data_out.build_patches ();
7637 * const std::string filename_base_stress = ("stress-polar-" + filename_base);
7639 * const std::string filename =
7640 * (output_dir + filename_base_stress + "-"
7641 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7643 * std::ofstream output_vtu((filename + ".vtu").c_str());
7644 * data_out.write_vtu(output_vtu);
7645 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7647 * if (this_mpi_process == 0)
7649 * std::vector<std::string> filenames;
7650 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7651 * filenames.push_back(filename_base_stress + "-" +
7652 * Utilities::int_to_string(i, 4) +
7655 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7656 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7658 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7659 * data_out.write_pvtu_record(visit_master_output, filenames);
7667 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7668 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7669 * in order to compute stresses on nodes (by applying nodal averaging)
7670 * Therefore, each vertex has one degree of freedom
7673 * FE_Q<dim> fe_1 (1);
7674 * DoFHandler<dim> dof_handler_1 (triangulation);
7675 * dof_handler_1.distribute_dofs (fe_1);
7677 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7678 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7680 * std::vector< std::vector< Vector<double> > >
7681 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7682 * for (unsigned int i=0; i<dim; ++i)
7683 * for (unsigned int j=0; j<dim; ++j)
7685 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7688 * Vector<double> counter_on_vertices (dof_handler_1.n_dofs());
7689 * counter_on_vertices = 0;
7691 * cell = dof_handler.begin_active();
7692 * dg_cell = history_dof_handler.begin_active();
7693 * typename DoFHandler<dim>::active_cell_iterator
7694 * cell_1 = dof_handler_1.begin_active();
7695 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7696 * if (cell->is_locally_owned())
7699 * for (unsigned int i=0; i<dim; ++i)
7700 * for (unsigned int j=0; j<dim; ++j)
7702 * dg_cell->get_dof_values (history_stress_field[i][j],
7703 * local_history_stress_fe_values[i][j]);
7706 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7708 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7713 * Point<dim> point1, point2;
7714 * point1 = cell_1->vertex(v);
7715 * point2 = dg_cell->vertex(v);
7716 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7723 * counter_on_vertices (dof_1_vertex) += 1;
7725 * for (unsigned int i=0; i<dim; ++i)
7726 * for (unsigned int j=0; j<dim; ++j)
7728 * history_stress_on_vertices[i][j](dof_1_vertex) +=
7729 * local_history_stress_fe_values[i][j](v);
7735 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7737 * for (unsigned int i=0; i<dim; ++i)
7738 * for (unsigned int j=0; j<dim; ++j)
7740 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7746 * DataOut<dim> data_out;
7747 * data_out.attach_dof_handler (dof_handler_1);
7750 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_rr_averaged");
7751 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_tt_averaged");
7752 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_rt_averaged");
7754 * data_out.build_patches ();
7756 * const std::string filename_base_stress = ("averaged-stress-polar-" + filename_base);
7758 * const std::string filename =
7759 * (output_dir + filename_base_stress + "-"
7760 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7762 * std::ofstream output_vtu((filename + ".vtu").c_str());
7763 * data_out.write_vtu(output_vtu);
7764 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7766 * if (this_mpi_process == 0)
7768 * std::vector<std::string> filenames;
7769 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7770 * filenames.push_back(filename_base_stress + "-" +
7771 * Utilities::int_to_string(i, 4) +
7774 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7775 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7777 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7778 * data_out.write_pvtu_record(visit_master_output, filenames);
7785 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7791 * if ( std::abs( (present_time/end_time)*(pressure/sigma_0) - 0.6 ) <
7792 * .501*(present_timestep/end_time)*(pressure/sigma_0) )
7797 * table_results_2: presenting the stress_rr and stress_tt on the nodes of bottom edge
7800 * const unsigned int face_id = 3;
7802 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
7804 * bool evaluation_face_found = false;
7806 * typename DoFHandler<dim>::active_cell_iterator
7807 * cell = dof_handler.begin_active(),
7808 * endc = dof_handler.end(),
7809 * cell_1 = dof_handler_1.begin_active();
7810 * for (; cell!=endc; ++cell, ++cell_1)
7811 * if (cell->is_locally_owned())
7813 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7815 * if (cell->face(face)->at_boundary()
7817 * cell->face(face)->boundary_id() == face_id)
7819 * if (!evaluation_face_found)
7821 * evaluation_face_found = true;
7825 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7827 * types::global_dof_index dof_1_vertex =
7828 * cell_1->face(face)->vertex_dof_index(v, 0);
7829 * if (!vertices_found[dof_1_vertex])
7832 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
7834 * table_results_2.add_value("x coordinate", vertex_coordinate[0]);
7835 * table_results_2.add_value("stress_rr", history_stress_on_vertices[0][0](dof_1_vertex));
7836 * table_results_2.add_value("stress_tt", history_stress_on_vertices[1][1](dof_1_vertex));
7837 * table_results_2.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7839 * vertices_found[dof_1_vertex] = true;
7848 * AssertThrow(evaluation_face_found, ExcInternalError());
7852 * table_results_3: presenting the mean stress_rr of the nodes on the inner radius
7855 * const unsigned int face_id_2 = 0;
7857 * Tensor<2, dim> stress_node,
7858 * mean_stress_polar;
7859 * mean_stress_polar = 0;
7861 * std::vector<bool> vertices_found_2 (dof_handler_1.n_dofs(), false);
7862 * unsigned int no_vertices_found = 0;
7864 * evaluation_face_found = false;
7866 * cell = dof_handler.begin_active(),
7867 * endc = dof_handler.end(),
7868 * cell_1 = dof_handler_1.begin_active();
7869 * for (; cell!=endc; ++cell, ++cell_1)
7870 * if (cell->is_locally_owned())
7872 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7874 * if (cell->face(face)->at_boundary()
7876 * cell->face(face)->boundary_id() == face_id_2)
7878 * if (!evaluation_face_found)
7880 * evaluation_face_found = true;
7884 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7886 * types::global_dof_index dof_1_vertex =
7887 * cell_1->face(face)->vertex_dof_index(v, 0);
7888 * if (!vertices_found_2[dof_1_vertex])
7890 * for (unsigned int ir=0; ir<dim; ++ir)
7891 * for (unsigned int ic=0; ic<dim; ++ic)
7892 * stress_node[ir][ic] = history_stress_on_vertices[ir][ic](dof_1_vertex);
7894 * mean_stress_polar += stress_node;
7896 * vertices_found_2[dof_1_vertex] = true;
7897 * ++no_vertices_found;
7906 * AssertThrow(evaluation_face_found, ExcInternalError());
7908 * mean_stress_polar /= no_vertices_found;
7910 * table_results_3.add_value("time step", timestep_no);
7911 * table_results_3.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7912 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
7913 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
7914 * table_results_3.add_value("radius", inner_radius);
7915 * table_results_3.add_value("mean stress_rr", mean_stress_polar[0][0]);
7916 * table_results_3.add_value("mean stress_tt", mean_stress_polar[1][1]);
7923 * else if (base_mesh == "Perforated_strip_tension")
7925 * const double imposed_displacement (0.00055),
7926 * inner_radius (0.05);
7931 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7932 * 3d and plane strain
7936 * table_results: Demonstrates the result of displacement at the top left corner versus imposed tension
7941 * const Point<dim> point_C(0., height);
7942 * Vector<double> disp_C(dim);
7946 * make a non-parallel copy of solution
7949 * Vector<double> copy_solution(solution);
7951 * typename Evaluation::PointValuesEvaluation<dim>::
7952 * PointValuesEvaluation point_values_evaluation(point_C);
7954 * point_values_evaluation.compute (dof_handler, copy_solution, disp_C);
7956 * table_results.add_value("time step", timestep_no);
7957 * table_results.add_value("Cells", triangulation.n_global_active_cells());
7958 * table_results.add_value("DoFs", dof_handler.n_dofs());
7959 * table_results.add_value("4*mu*u_C/(sigma_0*r)", 4*mu*disp_C(1)/(sigma_0*inner_radius));
7965 * compute average sigma_yy on the bottom edge
7968 * double stress_yy_av;
7971 * const unsigned int face_id = 1;
7973 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
7974 * unsigned int no_vertices_in_face = 0;
7976 * bool evaluation_face_found = false;
7978 * typename DoFHandler<dim>::active_cell_iterator
7979 * cell = dof_handler.begin_active(),
7980 * endc = dof_handler.end(),
7981 * cell_1 = dof_handler_1.begin_active();
7982 * for (; cell!=endc; ++cell, ++cell_1)
7983 * if (cell->is_locally_owned())
7985 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7987 * if (cell->face(face)->at_boundary()
7989 * cell->face(face)->boundary_id() == face_id)
7991 * if (!evaluation_face_found)
7993 * evaluation_face_found = true;
7997 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7999 * types::global_dof_index dof_1_vertex =
8000 * cell_1->face(face)->vertex_dof_index(v, 0);
8001 * if (!vertices_found[dof_1_vertex])
8003 * stress_yy_av += history_stress_on_vertices[1][1](dof_1_vertex);
8004 * ++no_vertices_in_face;
8006 * vertices_found[dof_1_vertex] = true;
8015 * AssertThrow(evaluation_face_found, ExcInternalError());
8017 * stress_yy_av /= no_vertices_in_face;
8023 * table_results_2: Demonstrate the stress_yy on the nodes of bottom edge
8027 * if ( std::abs( (stress_yy_av/sigma_0) - .91 ) < .2 )
8030 * if ( (timestep_no) % 19 == 0 )
8037 * const unsigned int face_id = 1;
8039 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
8041 * bool evaluation_face_found = false;
8043 * typename DoFHandler<dim>::active_cell_iterator
8044 * cell = dof_handler.begin_active(),
8045 * endc = dof_handler.end(),
8046 * cell_1 = dof_handler_1.begin_active();
8047 * for (; cell!=endc; ++cell, ++cell_1)
8048 * if (cell->is_locally_owned())
8050 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
8052 * if (cell->face(face)->at_boundary()
8054 * cell->face(face)->boundary_id() == face_id)
8056 * if (!evaluation_face_found)
8058 * evaluation_face_found = true;
8062 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
8064 * types::global_dof_index dof_1_vertex =
8065 * cell_1->face(face)->vertex_dof_index(v, 0);
8067 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
8069 * if (!vertices_found[dof_1_vertex] && std::abs(vertex_coordinate[2])<1.e-8)
8071 * table_results_2.add_value("x", vertex_coordinate[0]);
8072 * table_results_2.add_value("x/r", vertex_coordinate[0]/inner_radius);
8073 * table_results_2.add_value("stress_xx/sigma_0", history_stress_on_vertices[0][0](dof_1_vertex)/sigma_0);
8074 * table_results_2.add_value("stress_yy/sigma_0", history_stress_on_vertices[1][1](dof_1_vertex)/sigma_0);
8075 * table_results_2.add_value("stress_yy_av/sigma_0", stress_yy_av/sigma_0);
8076 * table_results_2.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8078 * vertices_found[dof_1_vertex] = true;
8087 * AssertThrow(evaluation_face_found, ExcInternalError());
8093 * table_results_3: Demonstrate the Stress_mean (average tensile stress)
8094 * on the bottom edge versus epsilon_yy on the bottom left corner
8098 * double strain_yy_A = 0.;
8102 * compute strain_yy_A
8103 * Since the point A is the node on the bottom left corner,
8104 * we need to work just with one element
8108 * const Point<dim> point_A(inner_radius, 0, 0);
8110 * Vector<double> local_strain_yy_values_at_qpoints (quadrature_formula.size()),
8111 * local_strain_yy_fe_values (history_fe.dofs_per_cell);
8113 * SymmetricTensor<2, dim> strain_at_qpoint;
8115 * typename DoFHandler<dim>::active_cell_iterator
8116 * cell = dof_handler.begin_active(),
8117 * endc = dof_handler.end(),
8118 * dg_cell = history_dof_handler.begin_active();
8120 * bool cell_found = false;
8122 * for (; cell!=endc; ++cell, ++dg_cell)
8123 * if (cell->is_locally_owned() && !cell_found)
8125 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
8126 * if ( std::fabs(cell->vertex(v)[0] - point_A[0])<1e-6 &&
8127 * std::fabs(cell->vertex(v)[1] - point_A[1])<1e-6 &&
8128 * std::fabs(cell->vertex(v)[2] - point_A[2])<1e-6)
8130 * PointHistory<dim> *local_quadrature_points_history
8131 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8132 * Assert (local_quadrature_points_history >=
8133 * &quadrature_point_history.front(),
8134 * ExcInternalError());
8135 * Assert (local_quadrature_points_history <
8136 * &quadrature_point_history.back(),
8137 * ExcInternalError());
8141 * Then loop over the quadrature points of this cell:
8144 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8146 * strain_at_qpoint = local_quadrature_points_history[q].old_strain;
8148 * local_strain_yy_values_at_qpoints(q) = strain_at_qpoint[1][1];
8151 * qpoint_to_dof_matrix.vmult (local_strain_yy_fe_values,
8152 * local_strain_yy_values_at_qpoints);
8154 * strain_yy_A = local_strain_yy_fe_values (v);
8156 * cell_found = true;
8164 * table_results_3.add_value("time step", timestep_no);
8165 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
8166 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
8167 * table_results_3.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8168 * table_results_3.add_value("mean_tensile_stress/sigma_0", stress_yy_av/sigma_0);
8169 * table_results_3.add_value("E*strain_yy-A/sigma_0", e_modulus*strain_yy_A/sigma_0);
8174 * if (std::abs(present_time-end_time) < 1.e-7)
8176 * table_results_2.set_precision("Imposed u_y", 6);
8177 * table_results_3.set_precision("Imposed u_y", 6);
8181 * else if (base_mesh == "Cantiliver_beam_3d")
8183 * const double pressure (6e6),
8189 * table_results: Demonstrates the result of displacement at the top front point, Point A
8195 * Quantity of interest:
8196 * displacement at Point A (x=0, y=height/2, z=length)
8202 * const Point<dim> point_A(0, height/2, length);
8203 * Vector<double> disp_A(dim);
8207 * make a non-parallel copy of solution
8210 * Vector<double> copy_solution(solution);
8212 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
8214 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
8216 * table_results.add_value("time step", timestep_no);
8217 * table_results.add_value("Cells", triangulation.n_global_active_cells());
8218 * table_results.add_value("DoFs", dof_handler.n_dofs());
8219 * table_results.add_value("pressure", pressure*present_time/end_time);
8220 * table_results.add_value("u_A", disp_A(1));
8226 * demonstrate the location and maximum von-Mises stress in the
8227 * specified domain close to the clamped face, z = 0
8228 * top domain: height/2 - thickness_flange <= y <= height/2
8229 * 0 <= z <= 2 * thickness_flange
8230 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
8231 * 0 <= z <= 2 * thickness_flange
8237 * double VM_stress_max (0);
8238 * Point<dim> point_max;
8240 * SymmetricTensor<2, dim> stress_at_qpoint;
8242 * typename DoFHandler<dim>::active_cell_iterator
8243 * cell = dof_handler.begin_active(),
8244 * endc = dof_handler.end();
8246 * const FEValuesExtractors::Vector displacement(0);
8248 * for (; cell!=endc; ++cell)
8249 * if (cell->is_locally_owned())
8251 * PointHistory<dim> *local_quadrature_points_history
8252 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8253 * Assert (local_quadrature_points_history >=
8254 * &quadrature_point_history.front(),
8255 * ExcInternalError());
8256 * Assert (local_quadrature_points_history <
8257 * &quadrature_point_history.back(),
8258 * ExcInternalError());
8262 * Then loop over the quadrature points of this cell:
8265 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8267 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
8269 * const double VM_stress = Evaluation::get_von_Mises_stress(stress_at_qpoint);
8270 * if (VM_stress > VM_stress_max)
8272 * VM_stress_max = VM_stress;
8273 * point_max = local_quadrature_points_history[q].point;
8279 * table_results.add_value("maximum von_Mises stress", VM_stress_max);
8280 * table_results.add_value("x", point_max[0]);
8281 * table_results.add_value("y", point_max[1]);
8282 * table_results.add_value("z", point_max[2]);
8295 * <a name="elastoplastic.cc-PlasticityContactProblemrun"></a>
8296 * <h4>PlasticityContactProblem::run</h4>
8300 * As in all other tutorial programs, the <code>run()</code> function contains
8301 * the overall logic. There is not very much to it here: in essence, it
8302 * performs the loops over all mesh refinement cycles, and within each, hands
8303 * things over to the Newton solver in <code>solve_newton()</code> on the
8304 * current mesh and calls the function that creates graphical output for
8305 * the so-computed solution. It then outputs some statistics concerning both
8306 * run times and memory consumption that has been collected over the course of
8307 * computations on this mesh.
8310 * template <int dim>
8312 * ElastoPlasticProblem<dim>::run ()
8314 * computing_timer.reset();
8317 * present_timestep = 1;
8325 * ----------------------------------------------------------------
8326 * base_mesh == "Thick_tube_internal_pressure"
8330 * const Point<dim> center(0, 0);
8331 * const double inner_radius = .1,
8332 * outer_radius = .2;
8334 * const SphericalManifold<dim> inner_boundary_description(center, inner_radius);
8335 * triangulation.set_manifold (0, inner_boundary_description);
8337 * const SphericalManifold<dim> outer_boundary_description(center, outer_radius);
8338 * triangulation.set_manifold (1, outer_boundary_description);
8342 * ----------------------------------------------------------------
8343 * base_mesh == "Perforated_strip_tension"
8347 * const double inner_radius = 0.05;
8349 * const CylinderBoundary<dim> inner_boundary_description(inner_radius, 2);
8350 * triangulation.set_manifold (10, inner_boundary_description);
8354 * ----------------------------------------------------------------
8360 * setup_quadrature_point_history ();
8362 * while (present_time < end_time)
8364 * present_time += present_timestep;
8367 * if (present_time > end_time)
8369 * present_timestep -= (present_time - end_time);
8370 * present_time = end_time;
8372 * pcout << std::endl;
8373 * pcout << "Time step " << timestep_no << " at time " << present_time
8376 * relative_error = max_relative_error * 10;
8377 * current_refinement_cycle = 0;
8384 * ------------------------ Refinement based on the relative error -------------------------------
8390 * while (relative_error >= max_relative_error)
8395 * if ( (timestep_no > 1) && (current_refinement_cycle>0) && (relative_error >= max_relative_error) )
8397 * pcout << "The relative error, " << relative_error
8398 * << " , is still more than maximum relative error, "
8399 * << max_relative_error << ", but we move to the next increment.\n";
8400 * relative_error = .1 * max_relative_error;
8403 * if (relative_error >= max_relative_error)
8405 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8406 * ++current_refinement_cycle;
8414 * ------------------------ Refinement based on the number of refinement --------------------------
8418 * bool continue_loop = true;
8419 * while (continue_loop)
8424 * if ( (timestep_no == 1) && (current_refinement_cycle < 1) )
8426 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8427 * ++current_refinement_cycle;
8431 * continue_loop = false;
8439 * -------------------------------------------------------------------------------------------------
8445 * solution += incremental_displacement;
8447 * update_quadrature_point_history ();
8449 * output_results((std::string("solution-") +
8450 * Utilities::int_to_string(timestep_no, 4)).c_str());
8452 * computing_timer.print_summary();
8453 * computing_timer.reset();
8455 * Utilities::System::MemoryStats stats;
8456 * Utilities::System::get_memory_stats(stats);
8457 * pcout << "Peak virtual memory used, resident in kB: " << stats.VmSize << " "
8458 * << stats.VmRSS << std::endl;
8461 * if (std::abs(present_time-end_time) < 1.e-7)
8463 * const std::string filename = (output_dir + "Results");
8465 * std::ofstream output_txt((filename + ".txt").c_str());
8467 * pcout << std::endl;
8468 * table_results.write_text(output_txt);
8469 * pcout << std::endl;
8470 * table_results_2.write_text(output_txt);
8471 * pcout << std::endl;
8472 * table_results_3.write_text(output_txt);
8473 * pcout << std::endl;
8478 * if (base_mesh == "Thick_tube_internal_pressure")
8480 * triangulation.reset_manifold (0);
8481 * triangulation.reset_manifold (1);
8483 * else if (base_mesh == "Perforated_strip_tension")
8485 * triangulation.reset_manifold (10);
8494 * <a name="elastoplastic.cc-Thecodemaincodefunction"></a>
8495 * <h3>The <code>main</code> function</h3>
8499 * There really isn't much to the <code>
main()</code> function. It looks
8503 *
int main (
int argc,
char *argv[])
8505 *
using namespace dealii;
8506 *
using namespace ElastoPlastic;
8512 *
const int dim = 3;
8513 *
ElastoPlasticProblem<dim>::declare_parameters(prm);
8516 *
std::cerr <<
"*** Call this program as <./elastoplastic input.prm>" << std::endl;
8520 *
prm.parse_input(argv[1]);
8523 *
ElastoPlasticProblem<dim> problem(prm);
8527 *
catch (std::exception &exc)
8529 *
std::cerr << std::endl << std::endl
8530 *
<<
"----------------------------------------------------"
8532 *
std::cerr <<
"Exception on processing: " << std::endl
8533 *
<< exc.what() << std::endl
8534 *
<<
"Aborting!" << std::endl
8535 *
<<
"----------------------------------------------------"
8542 *
std::cerr << std::endl << std::endl
8543 *
<<
"----------------------------------------------------"
8545 *
std::cerr <<
"Unknown exception!" << std::endl
8546 *
<<
"Aborting!" << std::endl
8547 *
<<
"----------------------------------------------------"
* * triangulation refine_global(n_initial_global_refinements)
* * for(const auto &cell :triangulation.active_cell_iterators())
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
* * const Number outer_radius
* * * struct InterferenceTaperTransform *
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
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
unsigned int depth_console(const unsigned int n)
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
void initialize(const SparsityPattern &sparsity_pattern)
numbers::NumberTraits< Number >::real_type norm() const
reference push_back(const T &value)
#define Assert(cond, exc)
#define DeclException1(Exception1, type1, outsequence)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
@ component_is_part_of_vector
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
constexpr types::blas_int zero
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
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 > 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)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
* * * ScaleZFunction< dim, Number, components >::ScaleZFunction * component(component)
* * * ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number > ThermoPlasticMaterial * mu(mu)
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, ShapeDataMemorySpace > shape_data, const ViewTypeIn in, ViewTypeOut out)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void load(Archive &ar, ::std_cxx26::inplace_vector< T, N > &vec, const unsigned int)
long double gamma(const unsigned int n)
int(&) functions(const void *v1, const void *v2)
constexpr types::boundary_id invalid_boundary_id
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
std_cxx26::inplace_vector< unsigned int, ReferenceCells::max_n_vertices< structdim >()> vertices
types::boundary_id boundary_id
std::vector< std::vector< bool > > constant_modes
double aggregation_threshold
unsigned int smoother_sweeps
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
1.9.8