159 * <a name=
"elastoplastic.cc-Includefiles"></a>
160 * <h3>Include files</h3>
161 * The
set of include files is not much of a surprise any more at
this time:
164 * #include <deal.II/base/conditional_ostream.h>
165 * #include <deal.II/base/parameter_handler.h>
166 * #include <deal.II/base/utilities.h>
167 * #include <deal.II/base/index_set.h>
168 * #include <deal.II/base/quadrature_lib.h>
169 * #include <deal.II/base/function.h>
170 * #include <deal.II/base/logstream.h>
171 * #include <deal.II/base/timer.h>
172 * #include <deal.II/base/table_handler.h>
174 * #include <deal.II/lac/vector.h>
175 * #include <deal.II/lac/full_matrix.h>
176 * #include <deal.II/lac/sparsity_tools.h>
177 * #include <deal.II/lac/sparse_matrix.h>
178 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
179 * #include <deal.II/lac/block_sparsity_pattern.h>
180 * #include <deal.II/lac/solver_bicgstab.h>
181 * #include <deal.II/lac/precondition.h>
182 * #include <deal.II/lac/affine_constraints.h>
183 * #include <deal.II/lac/trilinos_sparse_matrix.h>
184 * #include <deal.II/lac/trilinos_block_sparse_matrix.h>
185 * #include <deal.II/lac/trilinos_vector.h>
186 * #include <deal.II/lac/trilinos_precondition.h>
187 * #include <deal.II/lac/trilinos_solver.h>
188 * #include <deal.II/lac/sparse_direct.h>
190 * #include <deal.II/grid/tria.h>
191 * #include <deal.II/grid/grid_generator.h>
192 * #include <deal.II/grid/grid_refinement.h>
193 * #include <deal.II/grid/grid_tools.h>
194 * #include <deal.II/grid/tria_accessor.h>
195 * #include <deal.II/grid/tria_iterator.h>
196 * #include <deal.II/grid/grid_out.h>
197 * #include <deal.II/grid/manifold_lib.h>
199 * #include <deal.II/distributed/tria.h>
200 * #include <deal.II/distributed/grid_refinement.h>
201 * #include <deal.II/distributed/solution_transfer.h>
203 * #include <deal.II/dofs/dof_handler.h>
204 * #include <deal.II/dofs/dof_accessor.h>
205 * #include <deal.II/dofs/dof_renumbering.h>
206 * #include <deal.II/dofs/dof_tools.h>
208 * #include <deal.II/fe/fe_q.h>
209 * #include <deal.II/fe/fe_system.h>
210 * #include <deal.II/fe/fe_values.h>
211 * #include <deal.II/fe/fe_dgq.h>
212 * #include <deal.II/fe/fe_tools.h>
214 * #include <deal.II/numerics/vector_tools.h>
215 * #include <deal.II/numerics/matrix_tools.h>
216 * #include <deal.II/numerics/data_out.h>
217 * #include <deal.II/numerics/error_estimator.h>
218 * #include <deal.II/numerics/fe_field_function.h>
219 * #include <deal.II/numerics/solution_transfer.h>
223 * And here the only two
new things among the header files: an include file in
224 * which
symmetric tensors of rank 2 and 4 are implemented, as introduced in
228 * #include <deal.II/base/symmetric_tensor.h>
232 * And a header that implements filters
for iterators looping over all
233 * cells. We will use
this when selecting only those cells
for output that are
234 * owned by the present process in a %
parallel program:
237 * #include <deal.II/grid/filtered_iterator.h>
240 * #include <iostream>
244 * This
final include file provides the <code>mkdir</code> function
245 * that we will use to create a directory
for output files,
if necessary:
248 * #include <sys/stat.h>
250 *
namespace ElastoPlastic
256 *
const unsigned int n_slices,
257 *
const double height,
262 *
Assert (input.n_levels() == 1,
263 * ExcMessage (
"The input triangulations must be coarse meshes."));
266 *
Assert(result.n_cells()==0, ExcMessage(
"resultin Triangulation need to be empty upon calling extrude_triangulation."));
267 *
Assert(height>0, ExcMessage(
"The height in extrude_triangulation needs to be positive."));
268 *
Assert(n_slices>=2, ExcMessage(
"The number of slices in extrude_triangulation needs to be at least 2."));
270 * std::vector<Point<3> > points(n_slices*input.n_used_vertices());
271 * std::vector<CellData<3> > cells;
272 * cells.reserve((n_slices-1)*input.n_active_cells());
274 *
for (
unsigned int slice=0; slice<n_slices; ++slice)
276 *
for (
unsigned int i=0; i<input.n_vertices(); ++i)
279 *
if (input.get_used_vertices()[i])
281 *
const Point<2> &v = input.get_vertices()[i];
282 * points[i+slice*input.n_vertices()](0) = v(0);
283 * points[i+slice*input.n_vertices()](1) = v(1);
284 * points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1);
290 * cell = input.begin_active(); cell != input.end(); ++cell)
292 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
295 *
for (
unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
298 * = cell->vertex_index(v)+slice*input.n_used_vertices();
300 * = cell->vertex_index(v)+(slice+1)*input.n_used_vertices();
303 * this_cell.material_id = cell->material_id();
304 * cells.push_back(this_cell);
310 * s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
312 * cell = input.begin_active(); cell != input.end(); ++cell)
315 *
for (
unsigned int f=0; f<4; ++f)
316 *
if (cell->at_boundary(f))
319 * bid =
std::max(bid, quad.boundary_id);
320 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
322 * quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_used_vertices();
323 * quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_used_vertices();
324 * quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_used_vertices();
325 * quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_used_vertices();
326 * s.boundary_quads.push_back(quad);
332 * cell = input.begin_active(); cell != input.end(); ++cell)
336 * quad.vertices[0] = cell->vertex_index(0);
337 * quad.vertices[1] = cell->vertex_index(1);
338 * quad.vertices[2] = cell->vertex_index(2);
339 * quad.vertices[3] = cell->vertex_index(3);
340 * s.boundary_quads.push_back(quad);
342 * quad.boundary_id = bid + 2;
343 *
for (
int i=0; i<4; ++i)
344 * quad.vertices[i] += (n_slices-1)*input.n_used_vertices();
345 * s.boundary_quads.push_back(quad);
348 * result.create_triangulation (points,
353 *
namespace Evaluation
365 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
366 * + stress[1][1]*stress[1][1]
367 * - stress[0][0]*stress[1][1]
368 * + 3*stress[0][1]*stress[0][1]);
369 * }
else if (dim == 3)
371 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
372 * + stress[1][1]*stress[1][1]
373 * + stress[2][2]*stress[2][2]
374 * - stress[0][0]*stress[1][1]
375 * - stress[1][1]*stress[2][2]
376 * - stress[0][0]*stress[2][2]
377 * + 3*( stress[0][1]*stress[0][1]
378 * +stress[1][2]*stress[1][2]
379 * +stress[0][2]*stress[0][2]) );
384 * -----------------------------------------------
385 *
"Perforated_strip_tension"
387 *
const double von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
388 * + stress[1][1]*stress[1][1]
389 * - stress[0][0]*stress[1][1]
390 * + 3*stress[0][1]*stress[0][1]);
391 * -----------------------------------------------
393 * plane strain / 3
d case
399 * -----------------------------------------------
411 *
return von_Mises_stress;
416 *
class PointValuesEvaluation
419 * PointValuesEvaluation (
const Point<dim> &evaluation_point);
427 * <<
"The evaluation point " << arg1
428 * <<
" was not found among the vertices of the present grid.");
435 * PointValuesEvaluation<dim>::
436 * PointValuesEvaluation (
const Point<dim> &evaluation_point)
438 * evaluation_point (evaluation_point)
445 * PointValuesEvaluation<dim>::
450 *
const unsigned int dofs_per_vertex = dof_handler.get_fe().dofs_per_vertex;
452 * ExcDimensionMismatch (
point_values.size(), dofs_per_vertex));
456 * cell = dof_handler.begin_active(),
457 * endc = dof_handler.end();
458 *
bool evaluation_point_found =
false;
459 *
for (; (cell!=endc) && !evaluation_point_found; ++cell)
461 *
if (cell->is_locally_owned() && !evaluation_point_found)
462 *
for (
unsigned int vertex=0;
463 * vertex<GeometryInfo<dim>::vertices_per_cell;
466 *
if (cell->vertex(vertex).distance (evaluation_point)
468 * cell->diameter() * 1
e-8)
470 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
472 *
point_values[id] = solution(cell->vertex_dof_index(vertex,
id));
475 * evaluation_point_found =
true;
482 * ExcEvaluationPointNotFound(evaluation_point));
491 * <a name=
"elastoplastic.cc-ThecodePointHistorycodeclass"></a>
492 * <h3>The <code>PointHistory</code>
class</h3>
496 * As was mentioned in the introduction, we have to store the old stress in
497 * quadrature
point so that we can compute the residual forces at
this point
498 * during the next time step. This alone would not warrant a structure with
499 * only one member, but in more complicated applications, we would have to
500 * store more information in quadrature points as well, such as the history
501 * variables of plasticity, etc. In essence, we have to store everything
502 * that affects the present state of the material here, which in plasticity
503 * is determined by the deformation history variables.
507 * We will not give
this class any meaningful functionality beyond being
508 * able to store data, i.e. there are no constructors, destructors, or other
509 * member
functions. In such cases of `dumb
' classes, we usually opt to
510 * declare them as <code>struct</code> rather than <code>class</code>, to
511 * indicate that they are closer to C-style structures than C++-style
516 * struct PointHistory
518 * SymmetricTensor<2,dim> old_stress;
519 * SymmetricTensor<2,dim> old_strain;
527 * <a name="elastoplastic.cc-ThecodeConstitutiveLawcodeclasstemplate"></a>
528 * <h3>The <code>ConstitutiveLaw</code> class template</h3>
532 * This class provides an interface for a constitutive law, i.e., for the
533 * relationship between strain @f$\varepsilon(\mathbf u)@f$ and stress
534 * @f$\sigma@f$. In this example we are using an elastoplastic material behavior
535 * with linear, isotropic hardening. Such materials are characterized by
536 * Young's modulus @f$E@f$, Poisson
's ratio @f$\nu@f$, the initial yield stress
537 * @f$\sigma_0@f$ and the isotropic hardening parameter @f$\gamma@f$. For @f$\gamma =
538 * 0@f$ we obtain perfect elastoplastic behavior.
542 * As explained in the paper that describes this program, the first Newton
543 * steps are solved with a completely elastic material model to avoid having
544 * to deal with both nonlinearities (plasticity and contact) at once. To this
545 * end, this class has a function <code>set_sigma_0()</code> that we use later
546 * on to simply set @f$\sigma_0@f$ to a very large value -- essentially
547 * guaranteeing that the actual stress will not exceed it, and thereby
548 * producing an elastic material. When we are ready to use a plastic model, we
549 * set @f$\sigma_0@f$ back to its proper value, using the same function. As a
550 * result of this approach, we need to leave <code>sigma_0</code> as the only
551 * non-const member variable of this class.
555 * class ConstitutiveLaw
558 * ConstitutiveLaw (const double E,
560 * const double sigma_0,
561 * const double gamma);
564 * set_sigma_0 (double sigma_zero);
567 * get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
568 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
571 * get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
572 * const std::vector<Tensor<2, dim> > &point_hessian,
573 * Tensor<5, dim> &stress_strain_tensor_grad) const;
576 * get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor,
577 * SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
578 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
581 * const double kappa;
584 * const double gamma;
586 * const SymmetricTensor<4, dim> stress_strain_tensor_kappa;
587 * const SymmetricTensor<4, dim> stress_strain_tensor_mu;
592 * The constructor of the ConstitutiveLaw class sets the required material
593 * parameter for our deformable body. Material parameters for elastic
594 * isotropic media can be defined in a variety of ways, such as the pair @f$E,
595 * \nu@f$ (elastic modulus and Poisson's number),
using the Lame parameters
596 * @f$\lambda,mu@f$ or several other commonly used conventions. Here, the
597 * constructor takes a description of material parameters in the form of
598 * @f$E,\nu@f$, but since
this turns out to these are not the coefficients that
599 * appear in the equations of the plastic projector, we immediately convert
600 * them into the more suitable set @f$\kappa,\mu@f$ of bulk and shear moduli. In
601 * addition, the constructor takes @f$\sigma_0@f$ (the yield stress absent any
602 * plastic strain) and @f$\gamma@f$ (the hardening parameter) as arguments. In
603 *
this constructor, we also compute the two principal components of the
604 * stress-strain relation and its linearization.
608 * ConstitutiveLaw<dim>::ConstitutiveLaw (
double E,
615 * --------------------
617 * kappa (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (3 * (1 - 2 * (nu / (1+nu))))),
618 * mu (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))),
619 * --------------------
620 * 3
d and plane strain
623 * kappa (E / (3 * (1 - 2 * nu))),
624 * mu (E / (2 * (1 + nu))),
627 * --------------------
632 * stress_strain_tensor_kappa (kappa
635 * stress_strain_tensor_mu (2 * mu
644 * ConstitutiveLaw<dim>::set_sigma_0 (
double sigma_zero)
646 * sigma_0 = sigma_zero;
653 * <a name=
"elastoplastic.cc-ConstitutiveLawget_stress_strain_tensor"></a>
654 * <h4>ConstitutiveLaw::get_stress_strain_tensor</h4>
658 * This is the principal component of the constitutive law. It projects the
659 * deviatoric part of the stresses in a quadrature
point back to the yield
660 * stress (i.e., the original yield stress @f$\sigma_0@f$ plus the term that
661 * describes linear isotropic hardening). We need
this function to calculate
662 * the nonlinear residual in PlasticityContactProblem::residual_nl_system. The
663 * computations follow the formulas laid out in the introduction.
667 * The function returns whether the quadrature
point is plastic to allow
for
668 * some statistics
downstream on how many of the quadrature points are
669 * plastic and how many are elastic.
674 * ConstitutiveLaw<dim>::
679 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
685 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
688 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
690 * stress_strain_tensor = stress_strain_tensor_mu;
691 *
if (von_Mises_stress > sigma_0)
693 *
const double beta = sigma_0 / von_Mises_stress;
694 * stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
697 * stress_strain_tensor += stress_strain_tensor_kappa;
699 *
return (von_Mises_stress > sigma_0);
705 * ConstitutiveLaw<dim>::
711 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
714 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
716 *
if (von_Mises_stress > sigma_0)
719 *
const double deviator_strain_tensor_norm = deviator_strain_tensor.norm();
720 *
const double multiplier = -(1-
gamma)*sigma_0/(2*mu*
std::pow(deviator_strain_tensor_norm,3));
723 * multiplier_vector = 0;
725 *
for (
unsigned int i=0; i!=dim; ++i)
726 *
for (
unsigned int m=0; m!=dim; ++m)
727 *
for (
unsigned int n=0; n!=dim; ++n)
729 * multiplier_vector(i) += deviator_strain_tensor[m][n] *
730 * ( 0.5*( point_hessian[m][n][i] + point_hessian[n][m][i] )
731 * + ( m==n && dim==2 ? -1/dim*(point_hessian[0][0][i]
732 * + point_hessian[1][1][i]) : 0 )
733 * + ( m==n && dim==3 ? -1/dim*(point_hessian[0][0][i]
734 * + point_hessian[1][1][i]
735 * + point_hessian[2][2][i]) : 0 ) );
740 * -----------------------------------------------
741 *
"Perforated_strip_tension"
744 * -----------------------------------------------
746 * plane strain / 3
d case
749 *
const double VM_factor =
std::sqrt(1.5);
752 * -----------------------------------------------
758 *
for (
unsigned int i=0; i!=dim; ++i)
759 *
for (
unsigned int j=0; j!=dim; ++j)
760 *
for (
unsigned int k=0; k!=dim; ++k)
761 *
for (
unsigned int l=0;
l!=dim; ++
l)
762 *
for (
unsigned int m=0; m!=dim; ++m)
764 * stress_strain_tensor_grad[i][j][k][
l][m] = 1/VM_factor
766 * * stress_strain_tensor_mu[i][j][k][
l]
767 * * multiplier_vector(m);
773 * stress_strain_tensor_grad = 0;
776 *
return (von_Mises_stress > sigma_0);
783 * <a name=
"elastoplastic.cc-ConstitutiveLawget_linearized_stress_strain_tensors"></a>
784 * <h4>ConstitutiveLaw::get_linearized_stress_strain_tensors</h4>
788 * This function returns the linearized stress strain tensor, linearized
789 * around the solution @f$u^{i-1}@f$ of the previous Newton step @f$i-1@f$. The
790 * parameter <code>strain_tensor</code> (commonly denoted
791 * @f$\varepsilon(u^{i-1})@f$) must be passed as an argument, and serves as the
792 * linearization
point. The function returns the derivative of the nonlinear
793 * constitutive law in the variable stress_strain_tensor, as well as the
794 * stress-strain tensor of the linearized problem in
795 * stress_strain_tensor_linearized. See
796 * PlasticityContactProblem::assemble_nl_system where
this function is used.
801 * ConstitutiveLaw<dim>::
807 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
810 * stress_strain_tensor = stress_strain_tensor_mu;
811 * stress_strain_tensor_linearized = stress_strain_tensor_mu;
814 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
815 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
817 *
if (von_Mises_stress > sigma_0)
819 *
const double beta = sigma_0 / von_Mises_stress;
820 * stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
821 * stress_strain_tensor_linearized *= (
gamma + (1 -
gamma) * beta);
822 * deviator_stress_tensor /= deviator_stress_tensor_norm;
823 * stress_strain_tensor_linearized -= (1 -
gamma) * beta * 2 * mu
825 * deviator_stress_tensor);
828 * stress_strain_tensor += stress_strain_tensor_kappa;
829 * stress_strain_tensor_linearized += stress_strain_tensor_kappa;
834 * Finally, below we will need a function that computes the rotation
matrix
835 * induced by a displacement at a given
point. In fact, of course, the
836 * displacement at a single
point only has a direction and a magnitude, it
837 * is the change in direction and magnitude that induces rotations. In
839 * displacement, or, more specifically, from the curl.
843 * The formulas by which the rotation matrices are determined are a little
844 * awkward, especially in 3
d. For 2
d, there is a simpler way, so we
845 * implement
this function twice, once
for 2
d and once
for 3
d, so that we
846 * can compile and use the program in both space dimensions
if so desired --
847 * after all, deal.II is all about dimension
independent programming and
848 * reuse of algorithm thoroughly tested with cheap computations in 2
d,
for
849 * the more expensive computations in 3
d. Here is one
case, where we have to
850 * implement different algorithms
for 2
d and 3
d, but then can write the rest
851 * of the program in a way that is
independent of the space dimension.
855 * So, without further ado to the 2
d implementation:
859 * get_rotation_matrix (
const std::vector<
Tensor<1,2> > &grad_u)
863 * First, compute the curl of the velocity field from the
gradients. Note
864 * that we are in 2
d, so the rotation is a
scalar:
867 *
const double curl = (grad_u[1][0] - grad_u[0][1]);
871 * From
this, compute the
angle of rotation:
878 * And from
this, build the antisymmetric rotation
matrix:
881 *
const double t[2][2] = {{
cos(angle),
sin(angle) },
882 * {-
sin(angle),
cos(angle) }
890 * The 3
d case is a little more contrived:
894 * get_rotation_matrix (
const std::vector<
Tensor<1,3> > &grad_u)
898 * Again
first compute the curl of the velocity field. This time, it is a
902 *
const Point<3> curl (grad_u[2][1] - grad_u[1][2],
903 * grad_u[0][2] - grad_u[2][0],
904 * grad_u[1][0] - grad_u[0][1]);
908 * From
this vector,
using its magnitude, compute the tangent of the
angle
909 * of rotation, and from it the actual
angle:
912 *
const double tan_angle =
std::sqrt(curl*curl);
917 * Now, here
's one problem: if the angle of rotation is too small, that
918 * means that there is no rotation going on (for example a translational
919 * motion). In that case, the rotation matrix is the identity matrix.
923 * The reason why we stress that is that in this case we have that
924 * <code>tan_angle==0</code>. Further down, we need to divide by that
925 * number in the computation of the axis of rotation, and we would get
926 * into trouble when dividing doing so. Therefore, let's shortcut
this and
933 *
static const double rotation[3][3]
934 * = {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } };
941 * Otherwise compute the real rotation
matrix. The algorithm
for this is
942 * not exactly obvious, but can be found in a number of books,
943 * particularly on computer games where rotation is a very frequent
944 * operation. Online, you can find a description at
946 * the signs as here) at
952 *
const double t = 1-c;
954 *
const Point<3> axis = curl/tan_angle;
955 *
const double rotation[3][3]
957 * t *axis[0] *axis[0]+c,
958 * t *axis[0] *axis[1]+s *axis[2],
959 * t *axis[0] *axis[2]-s *axis[1]
962 * t *axis[0] *axis[1]-s *axis[2],
963 * t *axis[1] *axis[1]+c,
964 * t *axis[1] *axis[2]+s *axis[0]
967 * t *axis[0] *axis[2]+s *axis[1],
968 * t *axis[1] *axis[1]-s *axis[0],
969 * t *axis[2] *axis[2]+c
978 * <h3>Equation data: Body forces, boundary forces,
979 * incremental boundary values</h3>
983 * The following should be relatively standard. We need classes
for
984 * the boundary forcing term (which we here choose to be zero)
985 * and incremental boundary values.
988 *
namespace EquationData
1234 * ----------------------------- TimoshenkoBeam ---------------------------------------
1437 * ------------------------- Thick_tube_internal_pressure ----------------------------------
1606 * ------------------------- Perforated_strip_tension ----------------------------------
1776 * ------------------------- Cantiliver_beam_3d ----------------------------------
1779 *
template <
int dim>
1780 *
class IncrementalBoundaryForce :
public Function<dim>
1783 * IncrementalBoundaryForce (
const double present_time,
1784 *
const double end_time);
1796 *
const double present_time,
1802 *
template <
int dim>
1803 * IncrementalBoundaryForce<dim>::
1804 * IncrementalBoundaryForce (
const double present_time,
1805 *
const double end_time)
1808 * present_time (present_time),
1809 * end_time (end_time),
1814 *
template <
int dim>
1816 * IncrementalBoundaryForce<dim>::vector_value (
const Point<dim> &p,
1821 * ExcDimensionMismatch (values.size(), dim));
1823 *
const double eps = 1.e-7 * height;
1827 * pressure should be imposed on the top surface, y = height
1834 * values(1) = -pressure;
1838 * compute the fraction of imposed force
1841 *
const double frac = present_time/end_time;
1846 *
template <
int dim>
1848 * IncrementalBoundaryForce<dim>::
1849 * vector_value_list (
const std::vector<
Point<dim> > &points,
1852 *
const unsigned int n_points = points.size();
1854 *
Assert (value_list.size() == n_points,
1855 * ExcDimensionMismatch (value_list.size(), n_points));
1857 *
for (
unsigned int p=0; p<n_points; ++p)
1858 * IncrementalBoundaryForce<dim>::vector_value (points[p], value_list[p]);
1862 *
template <
int dim>
1866 * BodyForce () :
Functions::ZeroFunction<dim> (dim) {}
1870 *
template <
int dim>
1871 *
class IncrementalBoundaryValues :
public Function<dim>
1874 * IncrementalBoundaryValues (
const double present_time,
1875 *
const double end_time);
1888 *
const double present_time,
1893 *
template <
int dim>
1894 * IncrementalBoundaryValues<dim>::
1895 * IncrementalBoundaryValues (
const double present_time,
1896 *
const double end_time)
1899 * present_time (present_time),
1900 * end_time (end_time)
1904 *
template <
int dim>
1906 * IncrementalBoundaryValues<dim>::
1911 * ExcDimensionMismatch (values.size(), dim));
1918 *
template <
int dim>
1920 * IncrementalBoundaryValues<dim>::
1921 * vector_value_list (
const std::vector<
Point<dim> > &points,
1924 *
const unsigned int n_points = points.size();
1926 *
Assert (value_list.size() == n_points,
1927 * ExcDimensionMismatch (value_list.size(), n_points));
1929 *
for (
unsigned int p=0; p<n_points; ++p)
1930 * IncrementalBoundaryValues<dim>::vector_value (points[p], value_list[p]);
1935 * -------------------------------------------------------------------------------
1941 *
namespace DualFunctional
1944 *
template <
int dim>
1952 *
const ConstitutiveLaw<dim> &constitutive_law,
1958 *
template <
int dim>
1959 *
class PointValuesEvaluation :
public DualFunctionalBase<dim>
1962 * PointValuesEvaluation (
const Point<dim> &evaluation_point);
1968 *
const ConstitutiveLaw<dim> &constitutive_law,
1974 * <<
"The evaluation point " << arg1
1975 * <<
" was not found among the vertices of the present grid.");
1982 *
template <
int dim>
1983 * PointValuesEvaluation<dim>::
1984 * PointValuesEvaluation (
const Point<dim> &evaluation_point)
1986 * evaluation_point (evaluation_point)
1990 *
template <
int dim>
1992 * PointValuesEvaluation<dim>::
1995 *
const ConstitutiveLaw<dim> &,
1999 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2000 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2003 * cell_dual = dof_handler_dual.begin_active(),
2004 * endc_dual = dof_handler_dual.end();
2005 *
for (; cell_dual!=endc_dual; ++cell_dual)
2006 *
for (
unsigned int vertex=0;
2007 * vertex<GeometryInfo<dim>::vertices_per_cell;
2009 *
if (cell_dual->vertex(vertex).distance(evaluation_point)
2012 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
2014 * rhs_dual(cell_dual->vertex_dof_index(vertex,
id)) = 1;
2019 *
AssertThrow (
false, ExcEvaluationPointNotFound(evaluation_point));
2023 *
template <
int dim>
2024 *
class PointXDerivativesEvaluation :
public DualFunctionalBase<dim>
2027 * PointXDerivativesEvaluation (
const Point<dim> &evaluation_point);
2033 *
const ConstitutiveLaw<dim> &constitutive_law,
2039 * <<
"The evaluation point " << arg1
2040 * <<
" was not found among the vertices of the present grid.");
2047 *
template <
int dim>
2048 * PointXDerivativesEvaluation<dim>::
2049 * PointXDerivativesEvaluation (
const Point<dim> &evaluation_point)
2051 * evaluation_point (evaluation_point)
2055 *
template <
int dim>
2057 * PointXDerivativesEvaluation<dim>::
2060 *
const ConstitutiveLaw<dim> &,
2064 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2067 *
FEValues<dim> fe_values (dof_handler_dual.get_fe(), quadrature,
2071 *
const unsigned int n_q_points = fe_values.n_quadrature_points;
2072 *
Assert ( n_q_points==quadrature.size() , ExcInternalError() );
2073 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2076 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2078 *
double total_volume = 0;
2081 * cell = dof_handler_dual.begin_active(),
2082 * endc = dof_handler_dual.end();
2083 *
for (; cell!=endc; ++cell)
2084 *
if (cell->center().distance(evaluation_point) <=
2087 * fe_values.reinit (cell);
2090 *
for (
unsigned int q=0; q<n_q_points; ++q)
2092 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2094 * cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
2095 * fe_values.JxW (q);
2098 * total_volume += fe_values.JxW (q);
2101 * cell->get_dof_indices (local_dof_indices);
2102 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2104 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2109 * ExcEvaluationPointNotFound(evaluation_point));
2111 * rhs_dual *= 1./total_volume;
2116 *
template <
int dim>
2117 *
class MeanDisplacementFace :
public DualFunctionalBase<dim>
2120 * MeanDisplacementFace (
const unsigned int face_id,
2121 *
const std::vector<bool> comp_mask);
2127 *
const ConstitutiveLaw<dim> &constitutive_law,
2132 *
const unsigned int face_id;
2133 *
const std::vector<bool> comp_mask;
2137 *
template <
int dim>
2138 * MeanDisplacementFace<dim>::
2139 * MeanDisplacementFace (
const unsigned int face_id,
2140 *
const std::vector<bool> comp_mask )
2142 * face_id (face_id),
2143 * comp_mask (comp_mask)
2146 * ExcDimensionMismatch (comp_mask.size(), dim) );
2150 *
template <
int dim>
2152 * MeanDisplacementFace<dim>::
2155 *
const ConstitutiveLaw<dim> &,
2161 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2163 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2167 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2168 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2169 *
const unsigned int n_face_q_points = face_quadrature.size();
2172 * ExcDimensionMismatch (dofs_per_vertex, dim) );
2174 * std::vector<unsigned int> comp_vector(dofs_per_vertex);
2175 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2179 * comp_vector[i] = 1;
2185 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2189 * bound_size : size of the boundary, in 2
d is the length
2190 * and in the 3
d case, area
2193 *
double bound_size = 0.;
2196 * cell = dof_handler_dual.begin_active(),
2197 * endc = dof_handler_dual.end();
2198 *
bool evaluation_face_found =
false;
2199 *
for (; cell!=endc; ++cell)
2202 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2204 *
if (cell->face(face)->at_boundary()
2206 * cell->face(face)->boundary_id() == face_id)
2208 *
if (!evaluation_face_found)
2210 * evaluation_face_found =
true;
2212 * fe_face_values.reinit (cell, face);
2214 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2216 * bound_size += fe_face_values.JxW(q_point);
2218 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2220 *
const unsigned int
2221 * component_i = dof_handler_dual.get_fe().system_to_component_index(i).first;
2223 * cell_rhs(i) += (fe_face_values.shape_value(i,q_point) *
2224 * comp_vector[component_i] *
2225 * fe_face_values.JxW(q_point));
2233 * cell->get_dof_indices (local_dof_indices);
2234 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2236 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2241 *
AssertThrow(evaluation_face_found, ExcInternalError());
2243 * rhs_dual /= bound_size;
2248 *
template <
int dim>
2249 *
class MeanStressFace :
public DualFunctionalBase<dim>
2252 * MeanStressFace (
const unsigned int face_id,
2253 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2259 *
const ConstitutiveLaw<dim> &constitutive_law,
2264 *
const unsigned int face_id;
2265 *
const std::vector<std::vector<unsigned int> > comp_stress;
2269 *
template <
int dim>
2270 * MeanStressFace<dim>::
2271 * MeanStressFace (
const unsigned int face_id,
2272 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2274 * face_id (face_id),
2275 * comp_stress (comp_stress)
2278 * ExcDimensionMismatch (comp_stress.size(), dim) );
2282 *
template <
int dim>
2284 * MeanStressFace<dim>::
2287 *
const ConstitutiveLaw<dim> &constitutive_law,
2293 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2295 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2299 *
FEFaceValues<dim> fe_face_values_dual (dof_handler_dual.get_fe(), face_quadrature,
2302 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2303 *
const unsigned int n_face_q_points = face_quadrature.size();
2305 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_face_q_points);
2310 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2314 * bound_size : size of the boundary, in 2
d is the length
2315 * and in the 3
d case, area
2318 *
double bound_size = 0.;
2320 *
bool evaluation_face_found =
false;
2323 * cell_dual = dof_handler_dual.begin_active(),
2324 * endc_dual = dof_handler_dual.end(),
2325 * cell = dof_handler.begin_active();
2329 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2332 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2334 *
if (cell_dual->face(face)->at_boundary()
2336 * cell_dual->face(face)->boundary_id() == face_id)
2338 *
if (!evaluation_face_found)
2340 * evaluation_face_found =
true;
2343 * fe_face_values.reinit (cell, face);
2344 * fe_face_values_dual.reinit (cell_dual, face);
2346 * fe_face_values[displacement].get_function_symmetric_gradients(solution,
2349 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2351 * bound_size += fe_face_values_dual.JxW(q_point);
2353 * constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2354 * stress_strain_tensor);
2356 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2359 * stress_phi_i = stress_strain_tensor
2360 * * fe_face_values_dual[displacement].symmetric_gradient(i, q_point);
2362 *
for (
unsigned int k=0; k!=dim; ++k)
2364 *
for (
unsigned int l=0;
l!=dim; ++
l)
2366 *
if ( comp_stress[k][l] == 1 )
2368 * cell_rhs(i) += stress_phi_i[k][
l]
2370 * fe_face_values_dual.JxW(q_point);
2383 * cell_dual->get_dof_indices (local_dof_indices);
2384 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2386 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2391 *
AssertThrow(evaluation_face_found, ExcInternalError());
2393 * rhs_dual /= bound_size;
2398 *
template <
int dim>
2399 *
class MeanStressDomain :
public DualFunctionalBase<dim>
2402 * MeanStressDomain (
const std::string &base_mesh,
2403 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2409 *
const ConstitutiveLaw<dim> &constitutive_law,
2414 *
const std::string base_mesh;
2415 *
const std::vector<std::vector<unsigned int> > comp_stress;
2419 *
template <
int dim>
2420 * MeanStressDomain<dim>::
2421 * MeanStressDomain (
const std::string &base_mesh,
2422 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2424 * base_mesh (base_mesh),
2425 * comp_stress (comp_stress)
2428 * ExcDimensionMismatch (comp_stress.size(), dim) );
2432 *
template <
int dim>
2434 * MeanStressDomain<dim>::
2437 *
const ConstitutiveLaw<dim> &constitutive_law,
2441 *
AssertThrow (base_mesh ==
"Cantiliver_beam_3d", ExcNotImplemented());
2446 * Mean stress at the specified domain is of interest.
2447 * The interest domains are located on the bottom and top of the flanges
2448 * close to the clamped face, z = 0
2449 * top domain: height/2 - thickness_flange <= y <= height/2
2450 * 0 <= z <= 2 * thickness_flange
2451 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
2452 * 0 <= z <= 2 * thickness_flange
2458 *
const double height = 200
e-3,
2459 * thickness_flange = 10
e-3;
2461 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2463 *
const QGauss<dim> quadrature_formula(dof_handler_dual.get_fe().tensor_degree()+1);
2465 *
FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
2467 *
FEValues<dim> fe_values_dual (dof_handler_dual.get_fe(), quadrature_formula,
2470 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2471 *
const unsigned int n_q_points = quadrature_formula.size();
2473 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
2478 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2482 * domain_size : size of the interested domain, in 2
d is the area
2486 *
double domain_size = 0.;
2488 *
bool evaluation_domain_found =
false;
2491 * cell_dual = dof_handler_dual.begin_active(),
2492 * endc_dual = dof_handler_dual.end(),
2493 * cell = dof_handler.begin_active();
2497 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2499 *
const double y = cell->center()[1],
2500 * z = cell->center()[2];
2503 * top domain: height/2 - thickness_flange <= y <= height/2
2504 * 0 <= z <= 2 * thickness_flange
2505 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
2506 * 0 <= z <= 2 * thickness_flange
2509 *
if ( ((z > 0) && (z < 2*thickness_flange)) &&
2510 * ( ((y > height/2 - thickness_flange) && (y < height/2)) ||
2511 * ((y > -height/2) && (y < -height/2 + thickness_flange)) ) )
2515 *
if (!evaluation_domain_found)
2517 * evaluation_domain_found =
true;
2520 * fe_values.reinit(cell);
2521 * fe_values_dual.reinit(cell_dual);
2523 * fe_values[displacement].get_function_symmetric_gradients(solution,
2526 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
2528 * domain_size += fe_values_dual.JxW(q_point);
2530 * constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2531 * stress_strain_tensor);
2533 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2536 * stress_phi_i = stress_strain_tensor
2537 * * fe_values_dual[displacement].symmetric_gradient(i, q_point);
2539 *
for (
unsigned int k=0; k!=dim; ++k)
2541 *
for (
unsigned int l=0;
l!=dim; ++
l)
2543 *
if ( comp_stress[k][l] == 1 )
2545 * cell_rhs(i) += stress_phi_i[k][
l]
2547 * fe_values_dual.JxW(q_point);
2559 * cell_dual->get_dof_indices (local_dof_indices);
2560 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2562 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2567 *
AssertThrow(evaluation_domain_found, ExcInternalError());
2569 * rhs_dual /= domain_size;
2574 *
template <
int dim>
2575 *
class MeanStrainEnergyFace :
public DualFunctionalBase<dim>
2578 * MeanStrainEnergyFace (
const unsigned int face_id,
2582 *
void assemble_rhs_nonlinear (
const DoFHandler<dim> &primal_dof_handler,
2588 *
const unsigned int face_id;
2594 *
template <
int dim>
2595 * MeanStrainEnergyFace<dim>::
2596 * MeanStrainEnergyFace (
const unsigned int face_id,
2600 * face_id (face_id),
2601 * lambda_function (&lambda_function),
2602 * mu_function (&mu_function)
2606 *
template <
int dim>
2608 * MeanStrainEnergyFace<dim>::
2616 * Assemble right hand side of the dual problem when the quantity of interest is
2617 * a nonlinear functional. In
this case, the QoI should be linearized which depends
2618 * on the solution of the primal problem.
2619 * The extractor of the linearized QoI functional is the
gradient of the the original
2620 * QoI functional with the primal solution values.
2628 * rhs.reinit (dof_handler.n_dofs());
2630 *
const QGauss<dim-1> face_quadrature(dof_handler.get_fe().tensor_degree()+1);
2631 *
FEFaceValues<dim> primal_fe_face_values (primal_dof_handler.get_fe(), face_quadrature,
2639 *
const unsigned int dofs_per_vertex = primal_dof_handler.get_fe().dofs_per_vertex;
2640 *
const unsigned int n_face_q_points = face_quadrature.size();
2641 *
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
2644 * ExcDimensionMismatch (dofs_per_vertex, dim) );
2646 * std::vector< std::vector< Tensor<1,dim> > > primal_solution_gradients;
2647 * primal_solution_gradients.resize(n_face_q_points);
2649 * std::vector<std::vector<Tensor<2,dim> > > primal_solution_hessians;
2650 * primal_solution_hessians.resize (n_face_q_points);
2652 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
2654 * primal_solution_gradients[i].resize (dofs_per_vertex);
2655 * primal_solution_hessians[i].resize (dofs_per_vertex);
2658 * std::vector<double> lambda_values (n_face_q_points);
2659 * std::vector<double> mu_values (n_face_q_points);
2663 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2667 * bound_size : size of the boundary, in 2
d is the length
2668 * and in the 3
d case, area
2671 *
double bound_size = 0.;
2673 *
bool evaluation_face_found =
false;
2676 * primal_cell = primal_dof_handler.begin_active(),
2677 * primal_endc = primal_dof_handler.end();
2680 * cell = dof_handler.begin_active(),
2681 * endc = dof_handler.end();
2683 *
for (; cell!=endc; ++cell, ++primal_cell)
2686 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2688 *
if (cell->face(face)->at_boundary()
2690 * cell->face(face)->boundary_id() == face_id)
2692 *
if (!evaluation_face_found)
2694 * evaluation_face_found =
true;
2696 * primal_fe_face_values.reinit (primal_cell, face);
2698 * primal_fe_face_values.get_function_gradients (primal_solution,
2699 * primal_solution_gradients);
2701 * primal_fe_face_values.get_function_hessians (primal_solution,
2702 * primal_solution_hessians);
2704 * lambda_function->value_list (primal_fe_face_values.get_quadrature_points(), lambda_values);
2705 * mu_function->value_list (primal_fe_face_values.get_quadrature_points(), mu_values);
2707 * fe_face_values.reinit (cell, face);
2709 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2711 * bound_size += primal_fe_face_values.JxW(q_point);
2713 *
for (
unsigned int m=0; m<dofs_per_cell; ++m)
2715 *
const unsigned int
2716 * component_m = dof_handler.get_fe().system_to_component_index(m).first;
2718 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2720 *
for (
unsigned int j=0; j!=dofs_per_vertex; ++j)
2722 * cell_rhs(m) += fe_face_values.shape_value(m,q_point) *
2724 * lambda_values[q_point] *
2726 * primal_solution_hessians[q_point][i][i][component_m] * primal_solution_gradients[q_point][j][j]
2728 * primal_solution_gradients[q_point][i][i] * primal_solution_hessians[q_point][j][j][component_m]
2731 * mu_values[q_point] *
2733 * 2*primal_solution_hessians[q_point][j][i][component_m] * primal_solution_gradients[q_point][j][i]
2735 * primal_solution_hessians[q_point][i][j][component_m] * primal_solution_gradients[q_point][j][i]
2737 * primal_solution_gradients[q_point][i][j] * primal_solution_hessians[q_point][j][i][component_m]
2740 * primal_fe_face_values.JxW(q_point);
2753 * cell->get_dof_indices (local_dof_indices);
2754 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2756 * rhs(local_dof_indices[i]) += cell_rhs(i);
2761 *
AssertThrow(evaluation_face_found, ExcInternalError());
2763 * rhs *= 1./(2*bound_size);
2776 *
template <
int dim>
2783 *
const ConstitutiveLaw<dim> &constitutive_law,
2784 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2785 *
const unsigned int ×tep_no,
2786 *
const std::string &output_dir,
2787 *
const std::string &base_mesh,
2788 *
const double &present_time,
2789 *
const double &end_time);
2791 *
void compute_error_DWR (
Vector<float> &estimated_error_per_cell);
2796 *
void setup_system ();
2797 *
void compute_dirichlet_constraints ();
2798 *
void assemble_matrix ();
2799 *
void assemble_rhs ();
2801 *
void output_results ();
2807 *
const unsigned int fe_degree;
2810 *
const unsigned int fe_degree_dual;
2815 *
const QGauss<dim - 1> face_quadrature_formula;
2825 *
const ConstitutiveLaw<dim> constitutive_law;
2830 *
unsigned int timestep_no;
2831 * std::string output_dir;
2832 *
const std::string base_mesh;
2833 *
double present_time;
2843 *
const ConstitutiveLaw<dim> &constitutive_law,
2844 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2845 *
const unsigned int ×tep_no,
2846 *
const std::string &output_dir,
2847 *
const std::string &base_mesh,
2848 *
const double &present_time,
2849 *
const double &end_time)
2853 * solution(solution),
2854 * fe_degree(fe.tensor_degree()),
2855 * fe_degree_dual(fe_degree + 1),
2856 * fe_dual(
FE_Q<dim>(fe_degree_dual), dim),
2858 * quadrature_formula (fe_degree_dual + 1),
2859 * face_quadrature_formula (fe_degree_dual + 1),
2860 * constitutive_law (constitutive_law),
2862 * dual_functional (&dual_functional),
2863 * timestep_no (timestep_no),
2864 * output_dir (output_dir),
2865 * base_mesh (base_mesh),
2866 * present_time (present_time),
2867 * end_time (end_time)
2872 * DualSolver<dim>::~DualSolver()
2874 * dof_handler_dual.
clear ();
2879 *
void DualSolver<dim>::setup_system()
2881 * dof_handler.distribute_dofs(fe);
2883 * dof_handler_dual.distribute_dofs (fe_dual);
2884 * std::cout <<
" Number of degrees of freedom in dual problem: "
2885 * << dof_handler_dual.n_dofs()
2888 * constraints_hanging_nodes_dual.clear ();
2890 * constraints_hanging_nodes_dual);
2891 * constraints_hanging_nodes_dual.close ();
2893 * compute_dirichlet_constraints();
2895 * sparsity_pattern_dual.reinit (dof_handler_dual.n_dofs(),
2896 * dof_handler_dual.n_dofs(),
2897 * dof_handler_dual.max_couplings_between_dofs());
2902 * constraints_hanging_nodes_dual.condense (sparsity_pattern_dual);
2905 * constraints_dirichlet_and_hanging_nodes_dual.condense (sparsity_pattern_dual);
2907 * sparsity_pattern_dual.compress();
2909 * system_matrix_dual.reinit (sparsity_pattern_dual);
2911 * solution_dual.reinit (dof_handler_dual.n_dofs());
2912 * system_rhs_dual.reinit (dof_handler_dual.n_dofs());
2917 *
void DualSolver<dim>::compute_dirichlet_constraints()
2919 * constraints_dirichlet_and_hanging_nodes_dual.clear ();
2920 * constraints_dirichlet_and_hanging_nodes_dual.merge(constraints_hanging_nodes_dual);
2922 * std::vector<bool> component_mask(dim);
2924 *
if (base_mesh ==
"Timoshenko beam")
2928 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2929 * constraints_dirichlet_and_hanging_nodes_dual,
2932 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
2936 * the boundary x = 0
2939 * component_mask[0] =
true;
2940 * component_mask[1] =
false;
2943 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2944 * constraints_dirichlet_and_hanging_nodes_dual,
2948 * the boundary y = 0
2951 * component_mask[0] =
false;
2952 * component_mask[1] =
true;
2955 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2956 * constraints_dirichlet_and_hanging_nodes_dual,
2959 *
else if (base_mesh ==
"Perforated_strip_tension")
2963 * the boundary x = 0
2966 * component_mask[0] =
true;
2967 * component_mask[1] =
false;
2968 * component_mask[2] =
false;
2971 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2972 * constraints_dirichlet_and_hanging_nodes_dual,
2976 * the boundary y = 0
2979 * component_mask[0] =
false;
2980 * component_mask[1] =
true;
2981 * component_mask[2] =
false;
2984 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2985 * constraints_dirichlet_and_hanging_nodes_dual,
2989 * the boundary y = imposed incremental displacement
2992 * component_mask[0] =
false;
2993 * component_mask[1] =
true;
2994 * component_mask[2] =
false;
2997 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2998 * constraints_dirichlet_and_hanging_nodes_dual,
3001 *
else if (base_mesh ==
"Cantiliver_beam_3d")
3005 * the boundary x = y = z = 0
3008 * component_mask[0] =
true;
3009 * component_mask[1] =
true;
3010 * component_mask[2] =
true;
3013 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
3014 * constraints_dirichlet_and_hanging_nodes_dual,
3022 * constraints_dirichlet_and_hanging_nodes_dual.close();
3027 *
void DualSolver<dim>::assemble_matrix()
3034 *
const unsigned int dofs_per_cell_dual = fe_dual.dofs_per_cell;
3035 *
const unsigned int n_q_points = quadrature_formula.size();
3039 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell_dual);
3042 * cell_dual = dof_handler_dual.begin_active(),
3043 * endc_dual = dof_handler_dual.end(),
3044 * cell = dof_handler.begin_active();
3048 *
for (; cell_dual != endc_dual; ++cell_dual, ++cell)
3049 *
if (cell_dual->is_locally_owned())
3051 * fe_values.reinit(cell);
3053 * fe_values_dual.reinit(cell_dual);
3056 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3057 * fe_values[displacement].get_function_symmetric_gradients(solution,
3060 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3064 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3065 * stress_strain_tensor_linearized,
3066 * stress_strain_tensor);
3068 *
for (
unsigned int i = 0; i < dofs_per_cell_dual; ++i)
3071 * stress_phi_i = stress_strain_tensor_linearized
3072 * * fe_values_dual[displacement].symmetric_gradient(i, q_point);
3074 *
for (
unsigned int j = 0; j < dofs_per_cell_dual; ++j)
3076 * * fe_values_dual[displacement].symmetric_gradient(j, q_point)
3077 * * fe_values_dual.JxW(q_point));
3083 * cell_dual->get_dof_indices(local_dof_indices);
3084 * constraints_dirichlet_and_hanging_nodes_dual.distribute_local_to_global(cell_matrix,
3085 * local_dof_indices,
3086 * system_matrix_dual);
3094 *
void DualSolver<dim>::assemble_rhs()
3096 * dual_functional->assemble_rhs (dof_handler, solution, constitutive_law,
3097 * dof_handler_dual, system_rhs_dual);
3098 * constraints_dirichlet_and_hanging_nodes_dual.condense (system_rhs_dual);
3103 *
void DualSolver<dim>::solve()
3107 * +++ direct solver +++++++++
3115 * After the decomposition, we can use A_direct like a
matrix representing
3116 * the inverse of our system
matrix, so to compute the solution we just
3117 * have to multiply with the right hand side vector:
3120 * A_direct.vmult(solution_dual, system_rhs_dual);
3124 * ++++ iterative solver ++ CG ++++ doesn
't work
3125 * SolverControl solver_control (5000, 1e-12);
3126 * SolverCG<> cg (solver_control);
3130 * PreconditionSSOR<> preconditioner;
3131 * preconditioner.initialize(system_matrix_dual, 1.2);
3135 * cg.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3140 * ++++ iterative solver ++ BiCGStab ++++++ doesn't work
3147 * preconditioner.
initialize(system_matrix_dual, 1.0);
3151 * bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3156 * +++++++++++++++++++++++++++++++++++++++++++++++++
3162 * constraints_dirichlet_and_hanging_nodes_dual.distribute (solution_dual);
3166 *
void DualSolver<dim>::output_results()
3168 * std::string filename = (output_dir +
"dual-solution-" +
3170 * std::ofstream output (filename.c_str());
3173 * std::vector<std::string> solution_names;
3177 * solution_names.push_back (
"displacement");
3180 * solution_names.push_back (
"x_displacement");
3181 * solution_names.push_back (
"y_displacement");
3184 * solution_names.push_back (
"x_displacement");
3185 * solution_names.push_back (
"y_displacement");
3186 * solution_names.push_back (
"z_displacement");
3189 *
Assert (
false, ExcNotImplemented());
3191 * data_out.add_data_vector (solution_dual, solution_names);
3192 * data_out.build_patches ();
3193 * data_out.write_vtk (output);
3197 *
void DualSolver<dim>::compute_error_DWR (
Vector<float> &estimated_error_per_cell)
3204 * solve the dual problem
3208 * assemble_matrix ();
3211 * output_results ();
3215 * compuate the dual weights
3222 * constraints_dirichlet_and_hanging_nodes_dual,
3227 * constraints_hanging_nodes);
3228 * constraints_hanging_nodes.close();
3231 * constraints_dirichlet_and_hanging_nodes_dual,
3234 * constraints_hanging_nodes,
3239 * estimate the error
3249 *
const unsigned int n_q_points = quadrature_formula.size();
3250 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3254 * std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
3255 *
for (
unsigned int i=0; i!=n_q_points; ++i)
3257 * cell_hessians[i].resize (dim);
3259 * std::vector<Vector<double> > dual_weights_cell_values (n_q_points,
Vector<double>(dim));
3261 *
const EquationData::BodyForce<dim> body_force;
3262 * std::vector<Vector<double> > body_force_values (n_q_points,
Vector<double>(dim));
3272 * fe_face_values_neighbor (fe_dual, face_quadrature_formula,
3280 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
3281 * std::vector<Vector<double> > jump_residual (n_face_q_points,
Vector<double>(dim));
3282 * std::vector<Vector<double> > dual_weights_face_values (n_face_q_points,
Vector<double>(dim));
3284 * std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
3285 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3287 * cell_grads[i].resize (dim);
3289 * std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
3290 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3292 * neighbor_grads[i].resize (dim);
3300 *
typename std::map<typename DoFHandler<dim>::face_iterator,
Vector<double> >
3303 * cell = dof_handler_dual.begin_active(),
3304 * endc = dof_handler_dual.
end();
3305 *
for (; cell!=endc; ++cell)
3306 *
if (cell->is_locally_owned())
3308 *
for (
unsigned int face_no=0;
3309 * face_no<GeometryInfo<dim>::faces_per_cell;
3312 * face_integrals[cell->face(face_no)].reinit (dim);
3313 * face_integrals[cell->face(face_no)] = -1e20;
3317 * std::vector<Vector<float> > error_indicators_vector;
3323 * ----------------- estimate_some -------------------------
3326 * cell = dof_handler_dual.begin_active();
3327 *
unsigned int present_cell = 0;
3328 *
for (; cell!=endc; ++cell, ++present_cell)
3329 *
if (cell->is_locally_owned())
3333 * --------------- integrate_over_cell -------------------
3336 * fe_values.reinit(cell);
3337 * body_force.vector_value_list(fe_values.get_quadrature_points(),
3338 * body_force_values);
3339 * fe_values[displacement].get_function_symmetric_gradients(primal_solution,
3341 * fe_values.get_function_hessians(primal_solution, cell_hessians);
3343 * fe_values.get_function_values(dual_weights,
3344 * dual_weights_cell_values);
3346 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3348 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3349 * stress_strain_tensor_linearized,
3350 * stress_strain_tensor);
3351 * constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
3352 * cell_hessians[q_point],
3353 * stress_strain_tensor_grad);
3355 *
for (
unsigned int i=0; i!=dim; ++i)
3357 * error_indicators_vector[present_cell](i) +=
3358 * body_force_values[q_point](i)*
3359 * dual_weights_cell_values[q_point](i)*
3360 * fe_values.JxW(q_point);
3361 *
for (
unsigned int j=0; j!=dim; ++j)
3363 *
for (
unsigned int k=0; k!=dim; ++k)
3365 *
for (
unsigned int l=0;
l!=dim; ++
l)
3367 * error_indicators_vector[present_cell](i) +=
3368 * ( stress_strain_tensor[i][j][k][l]*
3369 * 0.5*(cell_hessians[q_point][k][l][j]
3371 * cell_hessians[q_point][l][k][j])
3372 * + stress_strain_tensor_grad[i][j][k][
l][j] * strain_tensor[q_point][k][
l]
3374 * dual_weights_cell_values[q_point](i) *
3375 * fe_values.JxW(q_point);
3385 * -------------------------------------------------------
3386 * compute face_integrals
3389 *
for (
unsigned int face_no=0;
3390 * face_no<GeometryInfo<dim>::faces_per_cell;
3393 *
if (cell->face(face_no)->at_boundary())
3395 *
for (
unsigned int id=0;
id!=dim; ++id)
3397 * face_integrals[cell->face(face_no)](id) = 0;
3402 *
if ((cell->neighbor(face_no)->has_children() ==
false) &&
3403 * (cell->neighbor(face_no)->level() == cell->level()) &&
3404 * (cell->neighbor(face_no)->index() < cell->index()))
3407 *
if (cell->at_boundary(face_no) ==
false)
3408 *
if (cell->neighbor(face_no)->level() < cell->level())
3412 *
if (cell->face(face_no)->has_children() ==
false)
3416 * ------------- integrate_over_regular_face -----------
3419 * fe_face_values_cell.reinit(cell, face_no);
3420 * fe_face_values_cell.get_function_gradients (primal_solution,
3424 * ExcInternalError());
3425 *
const unsigned int
3426 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3428 * neighbor = cell->neighbor(face_no);
3430 * fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
3431 * fe_face_values_neighbor.get_function_gradients (primal_solution,
3434 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3436 * q_cell_strain_tensor = 0.;
3437 * q_neighbor_strain_tensor = 0.;
3438 *
for (
unsigned int i=0; i!=dim; ++i)
3440 *
for (
unsigned int j=0; j!=dim; ++j)
3442 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3443 * cell_grads[q_point][j][i] );
3444 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3445 * neighbor_grads[q_point][j][i] );
3449 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3450 * cell_stress_strain_tensor);
3451 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3452 * neighbor_stress_strain_tensor);
3454 * jump_residual[q_point] = 0.;
3455 *
for (
unsigned int i=0; i!=dim; ++i)
3457 *
for (
unsigned int j=0; j!=dim; ++j)
3459 *
for (
unsigned int k=0; k!=dim; ++k)
3461 *
for (
unsigned int l=0;
l!=dim; ++
l)
3463 * jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
3464 * q_cell_strain_tensor[k][l]
3466 * neighbor_stress_strain_tensor[i][j][k][l]*
3467 * q_neighbor_strain_tensor[k][l] )*
3468 * fe_face_values_cell.normal_vector(q_point)[j];
3476 * fe_face_values_cell.get_function_values (dual_weights,
3477 * dual_weights_face_values);
3480 * face_integral_vector = 0;
3481 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3483 *
for (
unsigned int i=0; i!=dim; ++i)
3485 * face_integral_vector(i) += jump_residual[q_point](i) *
3486 * dual_weights_face_values[q_point](i) *
3487 * fe_face_values_cell.JxW(q_point);
3491 *
Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
3492 * ExcInternalError());
3494 *
for (
unsigned int i=0; i!=dim; ++i)
3496 *
Assert (face_integrals[cell->face(face_no)](i) == -1e20,
3497 * ExcInternalError());
3498 * face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
3504 * -----------------------------------------------------
3512 * ------------- integrate_over_irregular_face ---------
3516 * face = cell->face(face_no);
3518 * neighbor = cell->neighbor(face_no);
3520 * ExcInternalError());
3521 *
Assert (neighbor->has_children(),
3522 * ExcInternalError());
3524 *
const unsigned int
3525 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3527 *
for (
unsigned int subface_no=0;
3528 * subface_no<face->n_children(); ++subface_no)
3531 * neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
3532 *
Assert (neighbor_child->face(neighbor_neighbor) ==
3533 * cell->face(face_no)->child(subface_no),
3534 * ExcInternalError());
3536 * fe_subface_values_cell.reinit (cell, face_no, subface_no);
3537 * fe_subface_values_cell.get_function_gradients (primal_solution,
3539 * fe_face_values_neighbor.reinit (neighbor_child,
3540 * neighbor_neighbor);
3541 * fe_face_values_neighbor.get_function_gradients (primal_solution,
3544 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3546 * q_cell_strain_tensor = 0.;
3547 * q_neighbor_strain_tensor = 0.;
3548 *
for (
unsigned int i=0; i!=dim; ++i)
3550 *
for (
unsigned int j=0; j!=dim; ++j)
3552 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3553 * cell_grads[q_point][j][i] );
3554 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3555 * neighbor_grads[q_point][j][i] );
3559 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3560 * cell_stress_strain_tensor);
3561 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3562 * neighbor_stress_strain_tensor);
3564 * jump_residual[q_point] = 0.;
3565 *
for (
unsigned int i=0; i!=dim; ++i)
3567 *
for (
unsigned int j=0; j!=dim; ++j)
3569 *
for (
unsigned int k=0; k!=dim; ++k)
3571 *
for (
unsigned int l=0;
l!=dim; ++
l)
3573 * jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
3574 * q_cell_strain_tensor[k][l]
3576 * neighbor_stress_strain_tensor[i][j][k][l]*
3577 * q_neighbor_strain_tensor[k][l] )*
3578 * fe_face_values_neighbor.normal_vector(q_point)[j];
3586 * fe_face_values_neighbor.get_function_values (dual_weights,
3587 * dual_weights_face_values);
3590 * face_integral_vector = 0;
3591 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3593 *
for (
unsigned int i=0; i!=dim; ++i)
3595 * face_integral_vector(i) += jump_residual[q_point](i) *
3596 * dual_weights_face_values[q_point](i) *
3597 * fe_face_values_neighbor.JxW(q_point);
3601 *
for (
unsigned int i=0; i!=dim; ++i)
3603 * face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
3610 *
for (
unsigned int subface_no=0;
3611 * subface_no<face->n_children(); ++subface_no)
3613 *
Assert (face_integrals.find(face->child(subface_no)) !=
3614 * face_integrals.end(),
3615 * ExcInternalError());
3616 *
for (
unsigned int i=0; i!=dim; ++i)
3618 *
Assert (face_integrals[face->child(subface_no)](i) != -1e20,
3619 * ExcInternalError());
3620 *
sum(i) += face_integrals[face->child(subface_no)](i);
3623 *
for (
unsigned int i=0; i!=dim; ++i)
3625 * face_integrals[face](i) =
sum(i);
3631 * -----------------------------------------------------
3641 * ----------------------------------------------------------
3648 * cell = dof_handler_dual.begin_active();
3649 *
for (; cell!=endc; ++cell, ++present_cell)
3650 *
if (cell->is_locally_owned())
3652 *
for (
unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
3655 *
Assert(face_integrals.find(cell->face(face_no)) !=
3656 * face_integrals.end(),
3657 * ExcInternalError());
3659 *
for (
unsigned int id=0;
id!=dim; ++id)
3661 * error_indicators_vector[present_cell](id)
3662 * -= 0.5*face_integrals[cell->face(face_no)](id);
3667 * estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
3677 * <a name=
"elastoplastic.cc-ThecodePlasticityContactProblemcodeclasstemplate"></a>
3678 * <h3>The <code>PlasticityContactProblem</code>
class template</h3>
3682 * This is the main
class of this program and supplies all
functions
3683 * and variables needed to describe
3684 * the nonlinear contact problem. It is
3685 * close to @ref step_41
"step-41" but with some additional
3686 * features like handling hanging nodes,
3687 * a Newton method,
using Trilinos and p4est
3688 *
for parallel distributed computing.
3689 * To deal with hanging nodes makes
3690 * life a bit more complicated since
3692 * We create a Newton method
for the
3693 * active
set method
for the contact
3694 * situation and to handle the nonlinear
3695 *
operator for the constitutive law.
3699 * The
general layout of
this class is very much like for most other tutorial programs.
3700 * To make our life a bit easier,
this class reads a
set of input parameters from an input file. These
3701 * parameters,
using the
ParameterHandler class, are declared in the <code>declare_parameters</code>
3702 * function (which is
static so that it can be called before we even create an
object of the current
3703 * type), and a
ParameterHandler object that has been used to read an input file will then be passed
3704 * to the constructor of
this class.
3708 * The remaining member
functions are by and large as we have seen in several of the other tutorial
3709 * programs, though with additions
for the current nonlinear system. We will comment on their purpose
3710 * as we get to them further below.
3713 *
template <
int dim>
3714 *
class ElastoPlasticProblem
3724 *
void make_grid ();
3725 *
void setup_system ();
3726 *
void compute_dirichlet_constraints ();
3730 *
void solve_newton_system ();
3731 *
void solve_newton ();
3732 *
void compute_error ();
3734 *
void refine_grid ();
3736 *
void output_results (
const std::string &filename_base);
3740 * Next are three
functions that handle the history variables stored in each
3741 * quadrature
point. The
first one is called before the
first timestep to
3742 *
set up a pristine state
for the history variables. It only works on
3743 * those quadrature points on cells that belong to the present processor:
3746 *
void setup_quadrature_point_history ();
3750 * The
second one updates the history variables at the
end of each
3754 *
void update_quadrature_point_history ();
3758 * As far as member variables are concerned, we start with ones that we use to
3759 * indicate the
MPI universe
this program runs on, and then two
numbers
3760 * telling us how many participating processors there are, and where in
3761 *
this world we are., a stream we use to let
3762 * exactly one processor produce output to the console (see @ref step_17
"step-17") and
3763 * a variable that is used to time the various sections of the program:
3774 * The next
group describes the mesh and the finite element space.
3775 * In particular,
for this parallel program, the finite element
3776 * space has associated with it variables that indicate which degrees
3777 * of freedom live on the current processor (the index sets, see
3778 * also @ref step_40
"step-40" and the @ref distributed documentation module) as
3779 * well as a variety of constraints: those imposed by hanging nodes,
3780 * by Dirichlet boundary conditions, and by the active
set of
3782 * here, the
first only contains hanging node constraints, the
3783 *
second also those associated with Dirichlet boundary conditions,
3784 * and the third these plus the contact constraints.
3788 * The variable <code>active_set</code> consists of those degrees
3789 * of freedom constrained by the contact, and we use
3790 * <code>fraction_of_plastic_q_points_per_cell</code> to keep
3791 * track of the fraction of quadrature points on each cell where
3792 * the stress equals the yield stress. The latter is only used to
3793 * create graphical output showing the plastic zone, but not
for
3794 * any further computation; the variable is a member variable of
3795 *
this class since the information is computed as a by-product
3796 * of computing the residual, but is used only much later. (Note
3797 * that the vector is a vector of length
equal to the number of
3798 * active cells on the <i>local mesh</i>; it is never used to
3799 * exchange information between processors and can therefore be
3800 * a regular deal.II vector.)
3803 *
const unsigned int n_initial_global_refinements;
3806 *
const unsigned int fe_degree;
3820 * One difference of
this program is that we declare the quadrature
3821 * formula in the
class declaration. The reason is that in all the other
3822 * programs, it didn
't do much harm if we had used different quadrature
3823 * formulas when computing the matrix and the right hand side, for
3824 * example. However, in the present case it does: we store information in
3825 * the quadrature points, so we have to make sure all parts of the program
3826 * agree on where they are and how many there are on each cell. Thus, let
3827 * us first declare the quadrature formula that will be used throughout...
3830 * const QGauss<dim> quadrature_formula;
3831 * const QGauss<dim - 1> face_quadrature_formula;
3835 * ... and then also have a vector of history objects, one per quadrature
3836 * point on those cells for which we are responsible (i.e. we don't store
3837 * history data
for quadrature points on cells that are owned by other
3841 * std::vector<PointHistory<dim> > quadrature_point_history;
3845 * The way
this object is accessed is through a <code>user pointer</code>
3846 * that each cell, face, or edge holds: it is a <code>
void*</code> pointer
3847 * that can be used by application programs to associate arbitrary data to
3848 * cells, faces, or edges. What the program actually does with
this data
3849 * is within its own responsibility, the library just allocates some space
3850 *
for these pointers, and application programs can
set and read the
3851 * pointers
for each of these objects.
3858 * The next block of variables corresponds to the solution
3859 * and the linear systems we need to form. In particular,
this
3860 * includes the Newton
matrix and right hand side; the vector
3861 * that corresponds to the residual (i.e., the Newton right hand
3862 * side) but from which we have not eliminated the various
3863 * constraints and that is used to determine which degrees of
3864 * freedom need to be constrained in the next iteration; and
3865 * a vector that corresponds to the
diagonal of the @f$B@f$
matrix
3866 * briefly mentioned in the introduction and discussed in the
3867 * accompanying paper.
3879 * The next block of variables is then related to the time dependent
3880 * nature of the problem: they denote the length of the time interval
3881 * which we want to simulate, the present time and number of time step,
3882 * and length of present timestep:
3885 *
double present_time;
3886 *
double present_timestep;
3888 *
unsigned int timestep_no;
3892 * The next block contains the variables that describe the material
3896 *
const double e_modulus, nu, sigma_0,
gamma;
3897 * ConstitutiveLaw<dim> constitutive_law;
3901 * And then there is an assortment of other variables that are used
3902 * to identify the mesh we are asked to build as selected by the
3903 * parameter file, the obstacle that is being pushed into the
3904 * deformable body, the mesh refinement strategy, whether to transfer
3905 * the solution from one mesh to the next, and how many mesh
3906 * refinement cycles to perform. As possible, we mark these kinds
3907 * of variables as <code>
const</code> to help the reader identify
3908 * which ones may or may not be modified later on (the output directory
3909 * being an exception -- it is never modified outside the constructor
3910 * but it is awkward to initialize in the member-initializer-list
3911 * following the colon in the constructor since there we have only
3912 * one shot at setting it; the same is
true for the mesh refinement
3916 * const
std::string base_mesh;
3918 *
struct RefinementStrategy
3923 * refine_percentage,
3927 *
typename RefinementStrategy::value refinement_strategy;
3929 *
struct ErrorEstimationStrategy
3935 * weighted_residual_error,
3936 * weighted_kelly_error
3939 *
typename ErrorEstimationStrategy::value error_estimation_strategy;
3943 *
const bool transfer_solution;
3944 * std::string output_dir;
3949 *
unsigned int current_refinement_cycle;
3951 *
const double max_relative_error;
3952 *
float relative_error;
3954 *
const bool show_stresses;
3961 * <a name=
"elastoplastic.cc-ImplementationofthecodePlasticityContactProblemcodeclass"></a>
3962 * <h3>Implementation of the <code>PlasticityContactProblem</code>
class</h3>
3967 * <a name=
"elastoplastic.cc-PlasticityContactProblemdeclare_parameters"></a>
3968 * <h4>PlasticityContactProblem::declare_parameters</h4>
3972 * Let us start with the declaration of
run-time parameters that can be
3973 * selected in the input file. These values will be read back in the
3974 * constructor of
this class to initialize the member variables of this
3978 *
template <
int dim>
3982 * prm.declare_entry(
"polynomial degree",
"1",
3984 *
"Polynomial degree of the FE_Q finite element space, typically 1 or 2.");
3985 * prm.declare_entry(
"number of initial refinements",
"2",
3987 *
"Number of initial global mesh refinement steps before "
3988 *
"the first computation.");
3989 * prm.declare_entry(
"refinement strategy",
"percentage",
3991 *
"Mesh refinement strategy:\n"
3992 *
" global: one global refinement\n"
3993 *
" percentage: a fixed percentage of cells gets refined using the selected error estimator.");
3994 * prm.declare_entry(
"error estimation strategy",
"kelly_error",
3996 *
"Error estimation strategy:\n"
3997 *
" kelly_error: Kelly error estimator\n"
3998 *
" residual_error: residual-based error estimator\n"
3999 *
" weighted_residual_error: dual weighted residual (Goal-oriented) error estimator.\n");
4000 * prm.declare_entry(
"maximum relative error",
"0.05",
4002 *
"maximum relative error which plays the role of a criteria for refinement.");
4003 * prm.declare_entry(
"number of cycles",
"5",
4005 *
"Number of adaptive mesh refinement cycles to run.");
4006 * prm.declare_entry(
"output directory",
"",
4008 *
"Directory for output files (graphical output and benchmark "
4009 *
"statistics). If empty, use the current directory.");
4010 * prm.declare_entry(
"transfer solution",
"true",
4012 *
"Whether the solution should be used as a starting guess "
4013 *
"for the next finer mesh. If false, then the iteration starts at "
4014 *
"zero on every mesh.");
4015 * prm.declare_entry(
"base mesh",
"Thick_tube_internal_pressure",
4017 *
"Perforated_strip_tension|Cantiliver_beam_3d"),
4018 *
"Select the shape of the domain: 'box' or 'half sphere'");
4019 * prm.declare_entry(
"elasticity modulus",
"2.e11",
4021 *
"Elasticity modulus of the material in MPa (N/mm2)");
4022 * prm.declare_entry(
"Poissons ratio",
"0.3",
4024 *
"Poisson's ratio of the material");
4025 * prm.declare_entry(
"yield stress",
"2.e11",
4027 *
"Yield stress of the material in MPa (N/mm2)");
4028 * prm.declare_entry(
"isotropic hardening parameter",
"0.",
4030 *
"Isotropic hardening parameter of the material");
4031 * prm.declare_entry(
"show stresses",
"false",
4033 *
"Whether illustrates the stresses and von Mises stresses or not.");
4042 * <a name=
"elastoplastic.cc-ThecodePlasticityContactProblemcodeconstructor"></a>
4043 * <h4>The <code>PlasticityContactProblem</code> constructor</h4>
4047 * Given the declarations of member variables as well as the
4048 * declarations of
run-time parameters that are read from the input
4049 * file, there is
nothing surprising in
this constructor. In the body
4050 * we initialize the mesh refinement strategy and the output directory,
4051 * creating such a directory
if necessary.
4054 *
template <
int dim>
4055 * ElastoPlasticProblem<dim>::
4058 * mpi_communicator(MPI_COMM_WORLD),
4061 * pcout(std::cout, this_mpi_process == 0),
4065 * n_initial_global_refinements (prm.get_integer(
"number of initial refinements")),
4067 * fe_degree (prm.get_integer(
"polynomial degree")),
4070 * quadrature_formula (fe_degree + 1),
4071 * face_quadrature_formula (fe_degree + 1),
4073 * e_modulus (prm.get_double(
"elasticity modulus")),
4074 * nu (prm.get_double(
"Poissons ratio")),
4075 * sigma_0(prm.get_double(
"yield stress")),
4076 *
gamma (prm.get_double(
"isotropic hardening parameter")),
4077 * constitutive_law (e_modulus,
4082 * base_mesh (prm.get(
"base mesh")),
4084 * transfer_solution (prm.get_bool(
"transfer solution")),
4086 * table_results_2(),
4087 * table_results_3(),
4088 * max_relative_error (prm.get_double(
"maximum relative error")),
4089 * show_stresses (prm.get_bool(
"show stresses"))
4091 * std::string strat = prm.get(
"refinement strategy");
4092 *
if (strat ==
"global")
4093 * refinement_strategy = RefinementStrategy::refine_global;
4094 *
else if (strat ==
"percentage")
4095 * refinement_strategy = RefinementStrategy::refine_percentage;
4099 * strat = prm.get(
"error estimation strategy");
4100 *
if (strat ==
"kelly_error")
4101 * error_estimation_strategy = ErrorEstimationStrategy::kelly_error;
4102 *
else if (strat ==
"residual_error")
4103 * error_estimation_strategy = ErrorEstimationStrategy::residual_error;
4104 *
else if (strat ==
"weighted_residual_error")
4105 * error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error;
4109 * output_dir = prm.get(
"output directory");
4110 *
if (output_dir !=
"" && *(output_dir.rbegin()) !=
'/')
4111 * output_dir +=
"/";
4112 * mkdir(output_dir.c_str(), 0777);
4114 * pcout <<
" Using output directory '" << output_dir <<
"'" << std::endl;
4115 * pcout <<
" FE degree " << fe_degree << std::endl;
4116 * pcout <<
" transfer solution "
4117 * << (transfer_solution ?
"true" :
"false") << std::endl;
4125 * <a name=
"elastoplastic.cc-PlasticityContactProblemmake_grid"></a>
4126 * <h4>PlasticityContactProblem::make_grid</h4>
4130 * The next block deals with constructing the starting mesh.
4131 * We will use the following helper function and the
first
4132 * block of the <code>make_grid()</code> to construct a
4133 * mesh that corresponds to a half sphere. deal.II has a function
4134 * that creates such a mesh, but it is in the wrong location
4135 * and facing the wrong direction, so we need to
shift and
rotate
4136 * it a bit before
using it.
4140 * For later reference, as described in the documentation of
4142 * has boundary indicator zero,
while the remainder has boundary
4147 * rotate_half_sphere (
const Point<3> &in)
4149 *
return Point<3>(in(2), in(1), -in(0));
4152 *
template <
int dim>
4154 * ElastoPlasticProblem<dim>::make_grid ()
4156 *
if (base_mesh ==
"Timoshenko beam")
4160 *
const double length = .48,
4164 * point_2(length,
depth/2);
4166 * std::vector<unsigned int> repetitions(2);
4167 * repetitions[0] = 4;
4168 * repetitions[1] = 1;
4174 * give the indicators to boundaries
for specification,
4181 * |________________|
4183 * 0 to essential boundary conditions (left edge) which are as
default
4184 * 100 to the null boundaries (upper and lower edges) where we
do not need to take care of them
4185 * 5 to the natural boundaries (right edge)
for imposing the traction force
4191 *
for (; cell!=endc; ++cell)
4195 *
if ( std::fabs(cell->face(face)->center()(0)-length) < 1e-12 )
4197 * cell->face(face)->set_manifold_id(5);
4199 *
else if ( ( std::fabs(cell->face(face)->center()(1)-(
depth/2)) < 1e-12 )
4201 * ( std::fabs(cell->face(face)->center()(1)-(-
depth/2)) < 1e-12 ) )
4203 * cell->face(face)->set_manifold_id(100);
4212 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4216 * Example 1 from the paper: Zhong Z., .... A
new numerical method
for determining
4217 * collapse load-carrying capacity of structure made of elasto-plastic material,
4218 * J. Cent. South Univ. (2014) 21: 398-404
4224 *
const double inner_radius = .1,
4225 * outer_radius = .2;
4227 *
center, inner_radius, outer_radius,
4232 * give the indicators to boundaries
for specification,
4250 * 0 - inner boundary - natural boundary condition - impose the traction force
4251 * 1 - outer boundary -
free boundary - we
do not need to take care of them
4252 * 2 - left boundary - essential boundary condition - constrained to move along the x direction
4253 * 3 - bottom boundary - essential boundary condition - constrained to move along the y direction
4271 *
else if (base_mesh ==
"Perforated_strip_tension")
4275 * Example 2 from the paper: Zhong Z., .... A
new numerical method
for determining
4276 * collapse load-carrying capacity of structure made of elasto-plastic material,
4277 * J. Cent. South Univ. (2014) 21: 398-404
4282 *
const int dim_2d = 2;
4284 *
const double inner_radius = 0.05,
4285 * outer_radius = 0.1,
4287 * thickness = 0.004;
4300 *
const double eps = 1
e-7 * inner_radius;
4305 * center_2d, inner_radius, outer_radius,
4310 * Modify the triangulation_1
4314 * cell = triangulation_1.begin_active(),
4315 * endc = triangulation_1.end();
4316 * std::vector<bool> treated_vertices(triangulation_1.n_vertices(),
false);
4317 *
for (; cell != endc; ++cell)
4319 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4320 *
if (cell->face(f)->at_boundary() && cell->face(f)->center()(0)>eps &&
4321 * cell->face(f)->center()(1)>eps )
4328 *
point(0) = cell->face(f)->center()(0) - center_2d(0);
4329 *
point(1) = cell->face(f)->center()(1) - center_2d(1);
4332 *
for (
unsigned int v=0; v < GeometryInfo<dim_2d>::vertices_per_face; ++v)
4334 *
unsigned int vv = cell->face(f)->vertex_index(v);
4335 *
if (treated_vertices[vv] ==
false)
4337 * treated_vertices[vv] =
true;
4340 * cell->face(f)->vertex(v) = center_2d+
Point<dim_2d>(outer_radius,outer_radius);
4353 * Make the triangulation_2, a rectangular above the triangulation_1
4358 * point2 (outer_radius, height);
4366 * make the triangulation_2d and
refine it
4372 * Merge the two triangulation_1 and triangulation_2
4379 * Assign boundary indicators to the boundary faces
4400 * cell = triangulation_2d.begin_active(),
4401 * endc = triangulation_2d.end();
4402 *
for (; cell != endc; ++cell)
4404 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4406 *
if (cell->face(f)->at_boundary())
4408 *
if ( std::fabs(cell->face(f)->center()(1)) < eps )
4410 * cell->face(f)->set_manifold_id(1);
4412 *
else if ( std::fabs(cell->face(f)->center()(0)-outer_radius) < eps )
4414 * cell->face(f)->set_manifold_id(2);
4416 *
else if ( std::fabs(cell->face(f)->center()(1)-height) < eps )
4418 * cell->face(f)->set_manifold_id(3);
4420 *
else if ( std::fabs(cell->face(f)->center()(0)) < eps )
4422 * cell->face(f)->set_manifold_id(4);
4426 * cell->face(f)->set_all_boundary_ids(10);
4436 * triangulation_2d.set_manifold (10, inner_boundary_description);
4438 * triangulation_2d.refine_global(3);
4440 * triangulation_2d.reset_manifold (10);
4445 * Extrude the triangulation_2d and make it 3
d
4455 * Assign boundary indicators to the boundary faces
4481 *
for (; cell != endc; ++cell)
4483 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4485 *
if (cell->face(f)->at_boundary())
4487 * dist_vector = cell->face(f)->center() -
center;
4489 *
if ( std::fabs(dist_vector[1]) < eps )
4491 * cell->face(f)->set_manifold_id(1);
4493 *
else if ( std::fabs(dist_vector[0]-outer_radius) < eps )
4495 * cell->face(f)->set_manifold_id(2);
4497 *
else if ( std::fabs(dist_vector[1]-height) < eps )
4499 * cell->face(f)->set_manifold_id(3);
4501 *
else if ( std::fabs(dist_vector[0]) < eps )
4503 * cell->face(f)->set_manifold_id(4);
4505 *
else if ( std::fabs(dist_vector[2]) < eps )
4507 * cell->face(f)->set_manifold_id(5);
4509 *
else if ( std::fabs(dist_vector[2]-thickness) < eps )
4511 * cell->face(f)->set_manifold_id(6);
4515 * cell->face(f)->set_all_boundary_ids(10);
4532 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4536 * A rectangular tube made of Aluminium
4538 * approximation of beam 17250
4539 * units are in meter
4547 *
const int dim_2d = 2;
4549 *
const double length = .7,
4552 * thickness_web = 10
e-3,
4553 * thickness_flange = 10
e-3;
4561 *
const double eps = 1
e-7 * width;
4564 * Make the triangulation_b, a rectangular at the bottom of rectangular tube
4569 * point2 (width/2, -(height/2)+thickness_flange);
4571 * std::vector<unsigned int> repetitions(dim_2d);
4572 * repetitions[0] = 8;
4573 * repetitions[1] = 1;
4580 * Make the triangulation_t, a rectangular at the top of rectangular tube
4584 *
const Point<dim_2d> point1 (-width/2, (height/2)-thickness_flange),
4585 * point2 (width/2, height/2);
4587 * std::vector<unsigned int> repetitions(dim_2d);
4588 * repetitions[0] = 8;
4589 * repetitions[1] = 1;
4596 * Make the triangulation_l, a rectangular at the left of rectangular tube
4600 *
const Point<dim_2d> point1 (-width/2, -(height/2)+thickness_flange),
4601 * point2 (-(width/2)+thickness_web, (height/2)-thickness_flange);
4603 * std::vector<unsigned int> repetitions(dim_2d);
4604 * repetitions[0] = 1;
4605 * repetitions[1] = 18;
4612 * Make the triangulation_r, a rectangular at the right of rectangular tube
4616 *
const Point<dim_2d> point1 ((width/2)-thickness_web, -(height/2)+thickness_flange),
4617 * point2 (width/2, (height/2)-thickness_flange);
4619 * std::vector<unsigned int> repetitions(dim_2d);
4620 * repetitions[0] = 1;
4621 * repetitions[1] = 18;
4628 * make the triangulation_2d
4634 * merging every two triangles to make triangulation_2d
4638 * triangulation_blr;
4647 * Extrude the triangulation_2d and make it 3
d
4650 *
const unsigned int n_slices =
static_cast<int>(length*1000/20) + 1;
4656 * Assign boundary indicators to the boundary faces
4700 *
for (; cell != endc; ++cell)
4702 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4704 *
if (cell->face(f)->at_boundary())
4706 * dist_vector = cell->face(f)->center() -
center;
4708 *
if ( std::fabs(dist_vector[2]) < eps )
4710 * cell->face(f)->set_manifold_id(1);
4712 *
else if ( std::fabs(dist_vector[1]-(height/2)) < eps )
4714 * cell->face(f)->set_manifold_id(2);
4718 * cell->face(f)->set_all_boundary_ids(0);
4735 * pcout <<
" Number of active cells: "
4745 * <a name=
"elastoplastic.cc-PlasticityContactProblemsetup_system"></a>
4746 * <h4>PlasticityContactProblem::setup_system</h4>
4750 * The next piece in the puzzle is to
set up the
DoFHandler, resize
4751 * vectors and take care of various other status variables such as
4752 *
index sets and constraint matrices.
4756 * In the following, each
group of operations is put into a brace-enclosed
4757 * block that is being timed by the variable declared at the top of the
4759 * timed section, the destructor that is called at the end of the block
4763 *
template <
int dim>
4765 * ElastoPlasticProblem<dim>::setup_system ()
4771 * dof_handler.distribute_dofs(fe);
4772 * pcout <<
" Number of degrees of freedom: "
4773 * << dof_handler.n_dofs()
4776 * locally_owned_dofs = dof_handler.locally_owned_dofs();
4777 * locally_relevant_dofs =
4784 * constraints_hanging_nodes.reinit(locally_relevant_dofs);
4786 * constraints_hanging_nodes);
4787 * constraints_hanging_nodes.close();
4789 * pcout <<
" Number of active cells: "
4791 * <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
4794 * compute_dirichlet_constraints();
4800 *
if (timestep_no==1 || current_refinement_cycle!=0)
4802 * solution.reinit(locally_relevant_dofs, mpi_communicator);
4804 * incremental_displacement.reinit(locally_relevant_dofs, mpi_communicator);
4805 * newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
4806 * newton_rhs_residual.reinit(locally_owned_dofs, mpi_communicator);
4812 * Finally, we
set up sparsity patterns and matrices.
4813 * We temporarily (ab)use the system
matrix to also build the (diagonal)
4814 *
matrix that we use in eliminating degrees of freedom that are in contact
4815 * with the obstacle, but we then immediately
set the Newton
matrix back
4822 * mpi_communicator);
4825 * constraints_dirichlet_and_hanging_nodes,
false,
4826 * this_mpi_process);
4828 * newton_matrix.reinit(sp);
4836 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_dirichlet_constraints"></a>
4837 * <h4>PlasticityContactProblem::compute_dirichlet_constraints</h4>
4841 * This function, broken out of the preceding one, computes the constraints
4842 * associated with Dirichlet-type boundary conditions and puts them into the
4843 * <code>constraints_dirichlet_and_hanging_nodes</code> variable by merging
4844 * with the constraints that come from hanging nodes.
4848 * As laid out in the introduction, we need to distinguish between two
4850 * - If the domain is a box, we
set the displacement to zero at the bottom,
4851 * and allow vertical movement in z-direction along the sides. As
4852 * shown in the <code>make_grid()</code> function, the former corresponds
4853 * to boundary indicator 6, the latter to 8.
4854 * - If the domain is a half sphere, then we impose zero displacement along
4855 * the curved part of the boundary, associated with boundary indicator zero.
4858 *
template <
int dim>
4860 * ElastoPlasticProblem<dim>::compute_dirichlet_constraints ()
4862 * constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
4863 * constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
4865 * std::vector<bool> component_mask(dim);
4867 *
if (base_mesh ==
"Timoshenko beam")
4871 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4872 * constraints_dirichlet_and_hanging_nodes,
4875 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4879 * the boundary x = 0
4882 * component_mask[0] =
true;
4883 * component_mask[1] =
false;
4886 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4887 * constraints_dirichlet_and_hanging_nodes,
4891 * the boundary y = 0
4894 * component_mask[0] =
false;
4895 * component_mask[1] =
true;
4898 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4899 * constraints_dirichlet_and_hanging_nodes,
4902 *
else if (base_mesh ==
"Perforated_strip_tension")
4906 * the boundary x = 0
4909 * component_mask[0] =
true;
4910 * component_mask[1] =
false;
4911 * component_mask[2] =
false;
4914 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4915 * constraints_dirichlet_and_hanging_nodes,
4919 * the boundary y = 0
4922 * component_mask[0] =
false;
4923 * component_mask[1] =
true;
4924 * component_mask[2] =
false;
4927 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4928 * constraints_dirichlet_and_hanging_nodes,
4932 * the boundary y = imposed incremental displacement
4935 * component_mask[0] =
false;
4936 * component_mask[1] =
true;
4937 * component_mask[2] =
false;
4940 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4941 * constraints_dirichlet_and_hanging_nodes,
4944 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4948 * the boundary x = y = z = 0
4951 * component_mask[0] =
true;
4952 * component_mask[1] =
true;
4953 * component_mask[2] =
true;
4956 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4957 * constraints_dirichlet_and_hanging_nodes,
4966 * constraints_dirichlet_and_hanging_nodes.close();
4973 * <a name=
"elastoplastic.cc-PlasticityContactProblemassemble_newton_system"></a>
4974 * <h4>PlasticityContactProblem::assemble_newton_system</h4>
4978 * Given the complexity of the problem, it may come as a bit of a surprise
4979 * that assembling the linear system we have to solve in each Newton iteration
4980 * is actually fairly straightforward. The following function builds the Newton
4981 * right hand side and Newton
matrix. It looks fairly innocent because the
4982 * heavy lifting happens in the call to
4983 * <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
4985 * constraints we have previously computed.
4988 *
template <
int dim>
4990 * ElastoPlasticProblem<dim>::
4997 *
if (base_mesh ==
"Timoshenko beam")
4999 * traction_surface_id = 5;
5001 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5003 * traction_surface_id = 0;
5005 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5007 * traction_surface_id = 2;
5021 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5022 *
const unsigned int n_q_points = quadrature_formula.size();
5023 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5026 *
const EquationData::BodyForce<dim> body_force;
5027 * std::vector<Vector<double> > body_force_values(n_q_points,
5030 *
const EquationData::
5031 * IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5032 * std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5038 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5042 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5045 * std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
5048 * cell = dof_handler.begin_active(),
5049 * endc = dof_handler.end();
5053 *
for (; cell != endc; ++cell)
5054 *
if (cell->is_locally_owned())
5056 * fe_values.reinit(cell);
5060 * fe_values[displacement].get_function_symmetric_gradients(delta_linearization_point,
5061 * incremental_strain_tensor);
5065 * For assembling the local right hand side contributions, we need
5066 * to access the prior linearized stress
value in
this quadrature
5067 *
point. To get it, we use the user pointer of
this cell that
5068 * points into the global array to the quadrature
point data
5069 * corresponding to the
first quadrature
point of the present cell,
5070 * and then add an offset corresponding to the
index of the
5071 * quadrature
point we presently consider:
5074 *
const PointHistory<dim> *local_quadrature_points_history
5075 * =
reinterpret_cast<PointHistory<dim>*
>(cell->user_pointer());
5076 *
Assert (local_quadrature_points_history >=
5077 * &quadrature_point_history.front(),
5078 * ExcInternalError());
5079 *
Assert (local_quadrature_points_history <
5080 * &quadrature_point_history.back(),
5081 * ExcInternalError());
5085 * In addition, we need the values of the external body forces at
5086 * the quadrature points on
this cell:
5089 * body_force.vector_value_list(fe_values.get_quadrature_points(),
5090 * body_force_values);
5092 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5095 * tmp_strain_tensor_qpoint = local_quadrature_points_history[q_point].old_strain
5096 * + incremental_strain_tensor[q_point];
5100 * constitutive_law.get_linearized_stress_strain_tensors(tmp_strain_tensor_qpoint,
5101 * stress_strain_tensor_linearized,
5102 * stress_strain_tensor);
5105 *
for (
unsigned int i = 0; i < dim; ++i)
5107 * rhs_values_body_force[i] = body_force_values[q_point][i];
5110 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5114 * Having computed the stress-strain tensor and its linearization,
5115 * we can now put together the parts of the
matrix and right hand side.
5116 * In both, we need the linearized stress-strain tensor times the
5117 *
symmetric gradient of @f$\varphi_i@f$, i.e. the term @f$I_\Pi\varepsilon(\varphi_i)@f$,
5118 * so we introduce an abbreviation of
this term. Recall that the
5119 *
matrix corresponds to the bilinear form
5120 * @f$A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))@f$ in the
5121 * notation of the accompanying publication, whereas the right
5122 * hand side is @f$F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))@f$
5123 * where @f$u@f$ is the current linearization points (typically the last solution).
5124 * This might suggest that the right hand side will be zero
if the material
5125 * is completely elastic (where @f$I_\Pi=P_\Pi@f$) but
this ignores the fact
5126 * that the right hand side will also contain contributions from
5127 * non-homogeneous constraints due to the contact.
5131 * The code block that follows
this adds contributions that are due to
5132 * boundary forces, should there be any.
5136 * stress_phi_i = stress_strain_tensor_linearized
5137 * * fe_values[displacement].symmetric_gradient(i, q_point);
5139 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
5141 * * fe_values[displacement].symmetric_gradient(j, q_point)
5142 * * fe_values.JxW(q_point));
5146 * * incremental_strain_tensor[q_point] )
5148 * ( ( stress_strain_tensor
5149 * * fe_values[displacement].symmetric_gradient(i, q_point))
5150 * * tmp_strain_tensor_qpoint )
5152 * ( fe_values[displacement].
value(i, q_point)
5153 * * rhs_values_body_force )
5154 * ) * fe_values.JxW(q_point);
5159 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
5160 *
if (cell->face(face)->at_boundary()
5162 * cell->face(face)->boundary_id() == traction_surface_id)
5164 * fe_values_face.reinit(cell, face);
5166 * boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5167 * boundary_force_values);
5169 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
5172 *
for (
unsigned int i = 0; i < dim; ++i)
5174 * rhs_values[i] = boundary_force_values[q_point][i];
5176 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5177 * cell_rhs(i) += (fe_values_face[displacement].value(i, q_point)
5179 * * fe_values_face.JxW(q_point));
5183 * cell->get_dof_indices(local_dof_indices);
5184 * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_matrix, cell_rhs,
5185 * local_dof_indices,
5201 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_nonlinear_residual"></a>
5202 * <h4>PlasticityContactProblem::compute_nonlinear_residual</h4>
5206 * The following function computes the nonlinear residual of the equation
5207 * given the current solution (or any other linearization point). This
5208 * is needed in the linear search algorithm where we need to
try various
5209 * linear combinations of previous and current (trial) solution to
5210 * compute the (real, globalized) solution of the current Newton step.
5214 * That said, in a slight abuse of the name of the function, it actually
5215 * does significantly more. For example, it also computes the vector
5216 * that corresponds to the Newton residual but without eliminating
5217 * constrained degrees of freedom. We need
this vector to compute contact
5218 * forces and, ultimately, to compute the next active
set. Likewise, by
5219 * keeping track of how many quadrature points we encounter on each cell
5220 * that show plastic yielding, we also compute the
5221 * <code>fraction_of_plastic_q_points_per_cell</code> vector that we
5222 * can later output to visualize the plastic zone. In both of these cases,
5223 * the results are not necessary as part of the line search, and so we may
5224 * be wasting a small amount of time computing them. At the same time,
this
5225 * information appears as a natural by-product of what we need to
do here
5226 * anyway, and we want to collect it once at the
end of each Newton
5227 * step, so we may as well
do it here.
5231 * The actual implementation of
this function should be rather obvious:
5234 *
template <
int dim>
5236 * ElastoPlasticProblem<dim>::
5240 *
if (base_mesh ==
"Timoshenko beam")
5242 * traction_surface_id = 5;
5244 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5246 * traction_surface_id = 0;
5248 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5250 * traction_surface_id = 2;
5265 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5266 *
const unsigned int n_q_points = quadrature_formula.size();
5267 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5269 *
const EquationData::BodyForce<dim> body_force;
5270 * std::vector<Vector<double> > body_force_values(n_q_points,
5273 *
const EquationData::
5274 * IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5275 * std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5280 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5284 * newton_rhs_residual = 0;
5286 * fraction_of_plastic_q_points_per_cell = 0;
5289 * cell = dof_handler.begin_active(),
5290 * endc = dof_handler.end();
5291 *
unsigned int cell_number = 0;
5292 *
for (; cell != endc; ++cell, ++cell_number)
5293 *
if (cell->is_locally_owned())
5295 * fe_values.reinit(cell);
5298 * std::vector<SymmetricTensor<2, dim> > strain_tensors(n_q_points);
5299 * fe_values[displacement].get_function_symmetric_gradients(linearization_point,
5302 * body_force.vector_value_list(fe_values.get_quadrature_points(),
5303 * body_force_values);
5305 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5308 *
const bool q_point_is_plastic
5309 * = constitutive_law.get_stress_strain_tensor(strain_tensors[q_point],
5310 * stress_strain_tensor);
5311 *
if (q_point_is_plastic)
5312 * ++fraction_of_plastic_q_points_per_cell(cell_number);
5315 *
for (
unsigned int i = 0; i < dim; ++i)
5317 * rhs_values_body_force[i] = body_force_values[q_point][i];
5320 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5322 * cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5323 * * rhs_values_body_force
5325 * strain_tensors[q_point]
5326 * * stress_strain_tensor
5327 * * fe_values[displacement].symmetric_gradient(i, q_point)
5329 * * fe_values.JxW(q_point);
5333 * cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5335 * * fe_values.JxW(q_point));
5339 *
for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5340 *
if (cell->face(face)->at_boundary()
5341 * && cell->face(face)->boundary_id() == traction_surface_id)
5343 * fe_values_face.reinit(cell, face);
5345 * boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5346 * boundary_force_values);
5348 *
for (
unsigned int q_point = 0; q_point < n_face_q_points;
5352 *
for (
unsigned int i = 0; i < dim; ++i)
5354 * rhs_values[i] = boundary_force_values[q_point][i];
5356 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5357 * cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values
5358 * * fe_values_face.JxW(q_point));
5362 * cell->get_dof_indices(local_dof_indices);
5363 * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_rhs,
5364 * local_dof_indices,
5365 * newton_rhs_residual);
5369 * fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
5381 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton_system"></a>
5382 * <h4>PlasticityContactProblem::solve_newton_system</h4>
5386 * The last piece before we can discuss the actual Newton iteration
5387 * on a single mesh is the solver
for the linear systems. There are
5388 * a couple of complications that slightly obscure the code, but
5389 * mostly it is just setup then solve. Among the complications are:
5393 * - For the hanging nodes we have to apply
5395 * This is necessary
if a hanging node with solution
value @f$x_0@f$
5396 * has one neighbor with
value @f$x_1@f$ which is in contact with the
5397 * obstacle and one neighbor @f$x_2@f$ which is not in contact. Because
5398 * the update
for the former will be prescribed, the hanging node constraint
5399 * will have an inhomogeneity and will look like @f$x_0 = x_1/2 + \text{gap}/2@f$.
5400 * So the corresponding entries in the
5401 * ride-hang-side are non-zero with a
5402 * meaningless
value. These values we have to
5404 * - Like in @ref step_40
"step-40", we need to shuffle between vectors that
do and
do
5405 *
do not have ghost elements when solving or
using the solution.
5409 * The rest of the function is similar to @ref step_40
"step-40" and
5410 * @ref step_41
"step-41" except that we use a BiCGStab solver
5411 * instead of CG. This is due to the fact that
for very small hardening
5412 * parameters @f$\gamma@f$, the linear system becomes almost semidefinite though
5413 * still
symmetric. BiCGStab appears to have an easier time with such linear
5417 *
template <
int dim>
5419 * ElastoPlasticProblem<dim>::solve_newton_system ()
5424 * distributed_solution = incremental_displacement;
5426 * constraints_hanging_nodes.set_zero(distributed_solution);
5427 * constraints_hanging_nodes.set_zero(newton_rhs);
5431 * ------- Solver Bicgstab --- Preconditioner AMG -------------------
5438 * std::vector<std::vector<bool> > constant_modes;
5448 * additional_data.
w_cycle =
false;
5455 * preconditioner.initialize(newton_matrix, additional_data);
5470 *
const double relative_accuracy = 1
e-2;
5471 *
const double solver_tolerance = relative_accuracy
5472 * * newton_matrix.residual(tmp, distributed_solution,
5478 * solver_tolerance);
5480 * solver.solve(newton_matrix, distributed_solution,
5481 * newton_rhs, preconditioner);
5485 * pcout <<
" Error: " << solver_control.initial_value()
5486 * <<
" -> " << solver_control.last_value() <<
" in "
5487 * << solver_control.last_step() <<
" Bicgstab iterations."
5493 * ------- Solver CG --- Preconditioner SSOR -------------------
5501 * preconditioner.initialize(newton_matrix, additional_data);
5511 *
const double relative_accuracy = 1
e-8;
5514 *
const double relative_accuracy = 1
e-2;
5515 *
const double solver_tolerance = relative_accuracy
5516 * * newton_matrix.residual(tmp, distributed_solution,
5522 * solver_tolerance);
5526 * solver_tolerance);
5528 * solver.solve(newton_matrix, distributed_solution,
5529 * newton_rhs, preconditioner);
5531 * pcout <<
" Error: " << solver_control.initial_value()
5532 * <<
" -> " << solver_control.last_value() <<
" in "
5533 * << solver_control.last_step() <<
" CG iterations."
5538 * ........................................................
5544 * constraints_dirichlet_and_hanging_nodes.distribute(distributed_solution);
5546 * incremental_displacement = distributed_solution;
5553 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton"></a>
5554 * <h4>PlasticityContactProblem::solve_newton</h4>
5558 * This is,
finally, the function that implements the damped Newton method
5559 * on the current mesh. There are two nested loops: the outer
loop for the Newton
5560 * iteration and the inner
loop for the line search which
5561 * will be used only
if necessary. To obtain a good and reasonable
5562 * starting
value we solve an elastic problem in very
first Newton step on each
5563 * mesh (or only on the
first mesh
if we transfer solutions between meshes). We
5564 *
do so by setting the yield stress to an unreasonably large
value in these
5565 * iterations and then setting it back to the correct
value in subsequent
5570 * Other than
this, the top part of
this function should be reasonably
5574 *
template <
int dim>
5576 * ElastoPlasticProblem<dim>::solve_newton ()
5585 *
double residual_norm;
5586 *
double previous_residual_norm = -std::numeric_limits<double>::max();
5589 * previous_disp_norm = 0;
5591 *
const double correct_sigma = sigma_0;
5593 *
const unsigned int max_newton_iter = 100;
5595 *
for (
unsigned int newton_step = 1; newton_step <= max_newton_iter; ++newton_step)
5597 *
if (newton_step == 1
5599 * ((transfer_solution && timestep_no == 1)
5601 * !transfer_solution))
5602 * constitutive_law.set_sigma_0(1e+10);
5604 * constitutive_law.set_sigma_0(correct_sigma);
5606 * pcout <<
" " << std::endl;
5607 * pcout <<
" Newton iteration " << newton_step << std::endl;
5609 * pcout <<
" Assembling system... " << std::endl;
5610 * newton_matrix = 0;
5612 * newton_rhs_residual = 0;
5614 * tmp_solution = solution;
5615 * tmp_solution += incremental_displacement;
5616 * assemble_newton_system(tmp_solution,
5617 * incremental_displacement);
5619 * pcout <<
" Solving system... " << std::endl;
5620 * solve_newton_system();
5624 * It gets a bit more hairy after we have computed the
5625 * trial solution @f$\tilde{\mathbf u}@f$ of the current Newton step.
5626 * We handle a highly nonlinear problem so we have to damp
5627 * Newton
's method using a line search. To understand how we do this,
5628 * recall that in our formulation, we compute a trial solution
5629 * in each Newton step and not the update between old and new solution.
5630 * Since the solution set is a convex set, we will use a line
5631 * search that tries linear combinations of the
5632 * previous and the trial solution to guarantee that the
5633 * damped solution is in our solution set again.
5634 * At most we apply 5 damping steps.
5638 * There are exceptions to when we use a line search. First,
5639 * if this is the first Newton step on any mesh, then we don't have
5640 * any
point to compare the residual to, so we
always accept a full
5641 * step. Likewise,
if this is the
second Newton step on the
first mesh (or
5642 * the
second on any mesh
if we don
't transfer solutions from
5643 * mesh to mesh), then we have computed the first of these steps using
5644 * just an elastic model (see how we set the yield stress sigma to
5645 * an unreasonably large value above). In this case, the first Newton
5646 * solution was a purely elastic one, the second one a plastic one,
5647 * and any linear combination would not necessarily be expected to
5648 * lie in the feasible set -- so we just accept the solution we just
5653 * In either of these two cases, we bypass the line search and just
5654 * update residual and other vectors as necessary.
5657 * if ((newton_step==1)
5659 * (transfer_solution && newton_step == 2 && current_refinement_cycle == 0)
5661 * (!transfer_solution && newton_step == 2))
5663 * tmp_solution = solution;
5664 * tmp_solution += incremental_displacement;
5665 * compute_nonlinear_residual(tmp_solution);
5666 * old_solution = incremental_displacement;
5668 * residual = newton_rhs_residual;
5670 * residual.compress(VectorOperation::insert);
5672 * residual_norm = residual.l2_norm();
5674 * pcout << " Accepting Newton solution with residual: "
5675 * << residual_norm << std::endl;
5679 * for (unsigned int i = 0; i < 5; ++i)
5681 * distributed_solution = incremental_displacement;
5683 * const double alpha = std::pow(0.5, static_cast<double>(i));
5684 * tmp_vector = old_solution;
5685 * tmp_vector.sadd(1 - alpha, alpha, distributed_solution);
5687 * TimerOutput::Scope t(computing_timer, "Residual and lambda");
5689 * locally_relevant_tmp_vector = tmp_vector;
5690 * tmp_solution = solution;
5691 * tmp_solution += locally_relevant_tmp_vector;
5692 * compute_nonlinear_residual(tmp_solution);
5693 * residual = newton_rhs_residual;
5695 * residual.compress(VectorOperation::insert);
5697 * residual_norm = residual.l2_norm();
5699 * pcout << " Residual of the system: "
5700 * << residual_norm << std::endl
5701 * << " with a damping parameter alpha = " << alpha
5704 * if (residual_norm < previous_residual_norm)
5708 * incremental_displacement = tmp_vector;
5709 * old_solution = incremental_displacement;
5712 * disp_norm = incremental_displacement.l2_norm();
5717 * The final step is to check for convergence. If the residual is
5718 * less than a threshold of @f$10^{-10}@f$, then we terminate
5719 * the iteration on the current mesh:
5720 * if (residual_norm < 1e-10)
5723 * if (residual_norm < 1e-7)
5726 * pcout << " difference of two consecutive incremental displacement l2 norm : "
5727 * << std::abs(disp_norm - previous_disp_norm) << std::endl;
5728 * if ( std::abs(disp_norm - previous_disp_norm) < 1e-10 &&
5729 * (residual_norm < 1e-5 || std::abs(residual_norm - previous_residual_norm)<1e-9) )
5731 * pcout << " Convergence by difference of two consecutive solution! " << std::endl;
5736 * previous_residual_norm = residual_norm;
5737 * previous_disp_norm = disp_norm;
5744 * <a name="elastoplastic.cc-PlasticityContactProblemcompute_error"></a>
5745 * <h4>PlasticityContactProblem::compute_error</h4>
5751 * template <int dim>
5753 * ElastoPlasticProblem<dim>::compute_error ()
5755 * TrilinosWrappers::MPI::Vector tmp_solution(locally_owned_dofs, mpi_communicator);
5756 * tmp_solution = solution;
5757 * tmp_solution += incremental_displacement;
5759 * estimated_error_per_cell.reinit (triangulation.n_active_cells());
5760 * if (error_estimation_strategy == ErrorEstimationStrategy::kelly_error)
5762 * using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
5764 * KellyErrorEstimator<dim>::estimate(dof_handler,
5765 * QGauss<dim - 1>(fe.degree + 2),
5766 * std::map<types::boundary_id, const Function<dim> *>(),
5768 * estimated_error_per_cell);
5771 * else if (error_estimation_strategy == ErrorEstimationStrategy::residual_error)
5773 * compute_error_residual(tmp_solution);
5776 * else if (error_estimation_strategy == ErrorEstimationStrategy::weighted_residual_error)
5780 * make a non-parallel copy of tmp_solution
5783 * Vector<double> copy_solution(tmp_solution);
5787 * the dual function definition (it should be defined previously, e.g. input file)
5790 * if (base_mesh == "Timoshenko beam")
5792 * double length = .48,
5795 * const Point<dim> evaluation_point(length, -depth/2);
5797 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5799 * DualSolver<dim> dual_solver(triangulation, fe,
5801 * constitutive_law, dual_functional,
5802 * timestep_no, output_dir, base_mesh,
5803 * present_time, end_time);
5805 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5808 * else if (base_mesh == "Thick_tube_internal_pressure")
5810 * const unsigned int face_id = 0;
5811 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5812 * for (unsigned int i=0; i!=dim; ++i)
5814 * comp_stress[i].resize(dim);
5815 * for (unsigned int j=0; j!=dim; ++j)
5817 * comp_stress[i][j] = 1;
5821 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5823 * DualSolver<dim> dual_solver(triangulation, fe,
5825 * constitutive_law, dual_functional,
5826 * timestep_no, output_dir, base_mesh,
5827 * present_time, end_time);
5829 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5832 * else if (base_mesh == "Perforated_strip_tension")
5836 * .........................................
5837 * Mean stress_yy over the bottom boundary
5840 * const unsigned int face_id = 1;
5841 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5842 * for (unsigned int i=0; i!=dim; ++i)
5844 * comp_stress[i].resize(dim);
5845 * for (unsigned int j=0; j!=dim; ++j)
5847 * comp_stress[i][j] = 0;
5850 * comp_stress[1][1] = 1;
5852 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5856 * .........................................
5862 * DualSolver<dim> dual_solver(triangulation, fe,
5864 * constitutive_law, dual_functional,
5865 * timestep_no, output_dir, base_mesh,
5866 * present_time, end_time);
5868 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5871 * else if (base_mesh == "Cantiliver_beam_3d")
5875 * Quantity of interest:
5876 * -----------------------------------------------------------
5877 * displacement at Point A (x=0, y=height/2, z=length)
5881 * const double length = .7,
5884 * const Point<dim> evaluation_point(0, height/2, length);
5886 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5891 * -----------------------------------------------------------
5892 * Mean stress at the specified domain is of interest.
5893 * The interest domains are located on the bottom and top of the flanges
5894 * close to the clamped face, z = 0
5895 * top domain: height/2 - thickness_flange <= y <= height/2
5896 * 0 <= z <= 2 * thickness_flange
5897 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
5898 * 0 <= z <= 2 * thickness_flange
5904 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5905 * for (unsigned int i=0; i!=dim; ++i)
5907 * comp_stress[i].resize(dim);
5908 * for (unsigned int j=0; j!=dim; ++j)
5910 * comp_stress[i][j] = 1;
5913 * DualFunctional::MeanStressDomain<dim> dual_functional(base_mesh, comp_stress);
5917 * -----------------------------------------------------------
5923 * DualSolver<dim> dual_solver(triangulation, fe,
5925 * constitutive_law, dual_functional,
5926 * timestep_no, output_dir, base_mesh,
5927 * present_time, end_time);
5929 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5934 * AssertThrow(false, ExcNotImplemented());
5941 * AssertThrow(false, ExcNotImplemented());
5945 * relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm();
5947 * pcout << "Estimated relative error = " << relative_error << std::endl;
5951 * template <int dim>
5953 * ElastoPlasticProblem<dim>::compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution)
5955 * FEValues<dim> fe_values(fe, quadrature_formula,
5957 * update_gradients |
5959 * update_quadrature_points |
5960 * update_JxW_values);
5962 * const unsigned int n_q_points = quadrature_formula.size();
5963 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5964 * SymmetricTensor<4, dim> stress_strain_tensor_linearized;
5965 * SymmetricTensor<4, dim> stress_strain_tensor;
5966 * Tensor<5, dim> stress_strain_tensor_grad;
5967 * std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
5968 * for (unsigned int i=0; i!=n_q_points; ++i)
5970 * cell_hessians[i].resize (dim);
5972 * const EquationData::BodyForce<dim> body_force;
5974 * std::vector<Vector<double> > body_force_values (n_q_points, Vector<double>(dim));
5975 * const FEValuesExtractors::Vector displacement(0);
5978 * FEFaceValues<dim> fe_face_values_cell(fe, face_quadrature_formula,
5980 * update_quadrature_points|
5981 * update_gradients |
5982 * update_JxW_values |
5983 * update_normal_vectors),
5984 * fe_face_values_neighbor (fe, face_quadrature_formula,
5986 * update_gradients |
5987 * update_JxW_values |
5988 * update_normal_vectors);
5989 * FESubfaceValues<dim> fe_subface_values_cell (fe, face_quadrature_formula,
5990 * update_gradients);
5992 * const unsigned int n_face_q_points = face_quadrature_formula.size();
5993 * std::vector<Vector<double> > jump_residual (n_face_q_points, Vector<double>(dim));
5994 * std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
5995 * for (unsigned int i=0; i!=n_face_q_points; ++i)
5997 * cell_grads[i].resize (dim);
5999 * std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
6000 * for (unsigned int i=0; i!=n_face_q_points; ++i)
6002 * neighbor_grads[i].resize (dim);
6004 * SymmetricTensor<2, dim> q_cell_strain_tensor;
6005 * SymmetricTensor<2, dim> q_neighbor_strain_tensor;
6006 * SymmetricTensor<4, dim> cell_stress_strain_tensor;
6007 * SymmetricTensor<4, dim> neighbor_stress_strain_tensor;
6010 * typename std::map<typename DoFHandler<dim>::face_iterator, Vector<double> >
6012 * typename DoFHandler<dim>::active_cell_iterator
6013 * cell = dof_handler.begin_active(),
6014 * endc = dof_handler.end();
6015 * for (; cell!=endc; ++cell)
6016 * if (cell->is_locally_owned())
6018 * for (unsigned int face_no=0;
6019 * face_no<GeometryInfo<dim>::faces_per_cell;
6022 * face_integrals[cell->face(face_no)].reinit (dim);
6023 * face_integrals[cell->face(face_no)] = -1e20;
6027 * std::vector<Vector<float> > error_indicators_vector;
6028 * error_indicators_vector.resize( triangulation.n_active_cells(),
6029 * Vector<float>(dim) );
6033 * ----------------- estimate_some -------------------------
6036 * cell = dof_handler.begin_active();
6037 * unsigned int present_cell = 0;
6038 * for (; cell!=endc; ++cell, ++present_cell)
6039 * if (cell->is_locally_owned())
6043 * --------------- integrate_over_cell -------------------
6046 * fe_values.reinit(cell);
6047 * body_force.vector_value_list(fe_values.get_quadrature_points(),
6048 * body_force_values);
6049 * fe_values[displacement].get_function_symmetric_gradients(tmp_solution,
6051 * fe_values.get_function_hessians(tmp_solution, cell_hessians);
6053 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
6055 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
6056 * stress_strain_tensor_linearized,
6057 * stress_strain_tensor);
6058 * constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
6059 * cell_hessians[q_point],
6060 * stress_strain_tensor_grad);
6062 * for (unsigned int i=0; i!=dim; ++i)
6064 * error_indicators_vector[present_cell](i) +=
6065 * body_force_values[q_point](i)*fe_values.JxW(q_point);
6066 * for (unsigned int j=0; j!=dim; ++j)
6068 * for (unsigned int k=0; k!=dim; ++k)
6070 * for (unsigned int l=0; l!=dim; ++l)
6072 * error_indicators_vector[present_cell](i) +=
6073 * ( stress_strain_tensor[i][j][k][l]*
6074 * 0.5*(cell_hessians[q_point][k][l][j]
6076 * cell_hessians[q_point][l][k][j])
6077 * + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l]
6079 * fe_values.JxW(q_point);
6089 * -------------------------------------------------------
6090 * compute face_integrals
6093 * for (unsigned int face_no=0;
6094 * face_no<GeometryInfo<dim>::faces_per_cell;
6097 * if (cell->face(face_no)->at_boundary())
6099 * for (unsigned int id=0; id!=dim; ++id)
6101 * face_integrals[cell->face(face_no)](id) = 0;
6106 * if ((cell->neighbor(face_no)->has_children() == false) &&
6107 * (cell->neighbor(face_no)->level() == cell->level()) &&
6108 * (cell->neighbor(face_no)->index() < cell->index()))
6111 * if (cell->at_boundary(face_no) == false)
6112 * if (cell->neighbor(face_no)->level() < cell->level())
6116 * if (cell->face(face_no)->has_children() == false)
6120 * ------------- integrate_over_regular_face -----------
6123 * fe_face_values_cell.reinit(cell, face_no);
6124 * fe_face_values_cell.get_function_gradients (tmp_solution,
6127 * Assert (cell->neighbor(face_no).state() == IteratorState::valid,
6128 * ExcInternalError());
6129 * const unsigned int
6130 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6131 * const typename DoFHandler<dim>::active_cell_iterator
6132 * neighbor = cell->neighbor(face_no);
6134 * fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
6135 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6138 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6140 * q_cell_strain_tensor = 0.;
6141 * q_neighbor_strain_tensor = 0.;
6142 * for (unsigned int i=0; i!=dim; ++i)
6144 * for (unsigned int j=0; j!=dim; ++j)
6146 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6147 * cell_grads[q_point][j][i] );
6148 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6149 * neighbor_grads[q_point][j][i] );
6153 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6154 * cell_stress_strain_tensor);
6155 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6156 * neighbor_stress_strain_tensor);
6158 * jump_residual[q_point] = 0.;
6159 * for (unsigned int i=0; i!=dim; ++i)
6161 * for (unsigned int j=0; j!=dim; ++j)
6163 * for (unsigned int k=0; k!=dim; ++k)
6165 * for (unsigned int l=0; l!=dim; ++l)
6167 * jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
6168 * q_cell_strain_tensor[k][l]
6170 * neighbor_stress_strain_tensor[i][j][k][l]*
6171 * q_neighbor_strain_tensor[k][l] )*
6172 * fe_face_values_cell.normal_vector(q_point)[j];
6180 * Vector<double> face_integral_vector(dim);
6181 * face_integral_vector = 0;
6182 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6184 * for (unsigned int i=0; i!=dim; ++i)
6186 * face_integral_vector(i) += jump_residual[q_point](i) *
6187 * fe_face_values_cell.JxW(q_point);
6191 * Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
6192 * ExcInternalError());
6194 * for (unsigned int i=0; i!=dim; ++i)
6196 * Assert (face_integrals[cell->face(face_no)](i) == -1e20,
6197 * ExcInternalError());
6198 * face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
6204 * -----------------------------------------------------
6212 * ------------- integrate_over_irregular_face ---------
6215 * const typename DoFHandler<dim>::face_iterator
6216 * face = cell->face(face_no);
6217 * const typename DoFHandler<dim>::cell_iterator
6218 * neighbor = cell->neighbor(face_no);
6219 * Assert (neighbor.state() == IteratorState::valid,
6220 * ExcInternalError());
6221 * Assert (neighbor->has_children(),
6222 * ExcInternalError());
6224 * const unsigned int
6225 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6227 * for (unsigned int subface_no=0;
6228 * subface_no<face->n_children(); ++subface_no)
6230 * const typename DoFHandler<dim>::active_cell_iterator
6231 * neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
6232 * Assert (neighbor_child->face(neighbor_neighbor) ==
6233 * cell->face(face_no)->child(subface_no),
6234 * ExcInternalError());
6236 * fe_subface_values_cell.reinit (cell, face_no, subface_no);
6237 * fe_subface_values_cell.get_function_gradients (tmp_solution,
6239 * fe_face_values_neighbor.reinit (neighbor_child,
6240 * neighbor_neighbor);
6241 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6244 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6246 * q_cell_strain_tensor = 0.;
6247 * q_neighbor_strain_tensor = 0.;
6248 * for (unsigned int i=0; i!=dim; ++i)
6250 * for (unsigned int j=0; j!=dim; ++j)
6252 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6253 * cell_grads[q_point][j][i] );
6254 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6255 * neighbor_grads[q_point][j][i] );
6259 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6260 * cell_stress_strain_tensor);
6261 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6262 * neighbor_stress_strain_tensor);
6264 * jump_residual[q_point] = 0.;
6265 * for (unsigned int i=0; i!=dim; ++i)
6267 * for (unsigned int j=0; j!=dim; ++j)
6269 * for (unsigned int k=0; k!=dim; ++k)
6271 * for (unsigned int l=0; l!=dim; ++l)
6273 * jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
6274 * q_cell_strain_tensor[k][l]
6276 * neighbor_stress_strain_tensor[i][j][k][l]*
6277 * q_neighbor_strain_tensor[k][l] )*
6278 * fe_face_values_neighbor.normal_vector(q_point)[j];
6286 * Vector<double> face_integral_vector(dim);
6287 * face_integral_vector = 0;
6288 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6290 * for (unsigned int i=0; i!=dim; ++i)
6292 * face_integral_vector(i) += jump_residual[q_point](i) *
6293 * fe_face_values_neighbor.JxW(q_point);
6297 * for (unsigned int i=0; i!=dim; ++i)
6299 * face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
6304 * Vector<double> sum (dim);
6306 * for (unsigned int subface_no=0;
6307 * subface_no<face->n_children(); ++subface_no)
6309 * Assert (face_integrals.find(face->child(subface_no)) !=
6310 * face_integrals.end(),
6311 * ExcInternalError());
6312 * for (unsigned int i=0; i!=dim; ++i)
6314 * Assert (face_integrals[face->child(subface_no)](i) != -1e20,
6315 * ExcInternalError());
6316 * sum(i) += face_integrals[face->child(subface_no)](i);
6319 * for (unsigned int i=0; i!=dim; ++i)
6321 * face_integrals[face](i) = sum(i);
6327 * -----------------------------------------------------
6337 * ----------------------------------------------------------
6344 * cell = dof_handler.begin_active();
6345 * for (; cell!=endc; ++cell, ++present_cell)
6346 * if (cell->is_locally_owned())
6348 * for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
6351 * Assert(face_integrals.find(cell->face(face_no)) !=
6352 * face_integrals.end(),
6353 * ExcInternalError());
6355 * for (unsigned int id=0; id!=dim; ++id)
6357 * error_indicators_vector[present_cell](id)
6358 * -= 0.5*face_integrals[cell->face(face_no)](id);
6363 * estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
6373 * <a name="elastoplastic.cc-PlasticityContactProblemrefine_grid"></a>
6374 * <h4>PlasticityContactProblem::refine_grid</h4>
6378 * If you've made it
this far into the deal.II tutorial, the following
6379 * function refining the mesh should not pose any challenges to you
6380 * any more. It refines the mesh, either globally or
using the Kelly
6381 * error estimator, and
if so asked also transfers the solution from
6382 * the previous to the next mesh. In the latter
case, we also need
6383 * to compute the active set and other quantities again,
for which we
6384 * need the information computed by <code>compute_nonlinear_residual()</code>.
6387 * template <int dim>
6389 * ElastoPlasticProblem<dim>::refine_grid ()
6393 * ---------------------------------------------------------------
6394 * Make a field variable
for history variables to be able to
6395 * transfer the data to the quadrature points of the
new mesh
6400 * history_dof_handler.distribute_dofs (history_fe);
6401 * std::vector< std::vector< Vector<double> > >
6402 * history_stress_field (dim, std::vector<
Vector<double> >(dim)),
6403 * local_history_stress_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6404 * local_history_stress_fe_values (dim, std::vector<
Vector<double> >(dim));
6407 * std::vector< std::vector< Vector<double> > >
6408 * history_strain_field (dim, std::vector<
Vector<double> >(dim)),
6409 * local_history_strain_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6410 * local_history_strain_fe_values (dim, std::vector<
Vector<double> >(dim));
6412 *
for (
unsigned int i=0; i<dim; ++i)
6413 *
for (
unsigned int j=0; j<dim; ++j)
6415 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6416 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6417 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6419 * history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6420 * local_history_strain_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6421 * local_history_strain_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6424 * quadrature_formula.size());
6427 * quadrature_formula, quadrature_formula,
6428 * qpoint_to_dof_matrix);
6430 * cell = dof_handler.begin_active(),
6431 * endc = dof_handler.end(),
6432 * dg_cell = history_dof_handler.begin_active();
6433 *
for (; cell!=endc; ++cell, ++dg_cell)
6434 *
if (cell->is_locally_owned())
6436 * PointHistory<dim> *local_quadrature_points_history
6437 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6438 *
Assert (local_quadrature_points_history >=
6439 * &quadrature_point_history.front(),
6440 * ExcInternalError());
6441 *
Assert (local_quadrature_points_history <
6442 * &quadrature_point_history.back(),
6443 * ExcInternalError());
6444 *
for (
unsigned int i=0; i<dim; ++i)
6445 *
for (
unsigned int j=0; j<dim; ++j)
6447 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6449 * local_history_stress_values_at_qpoints[i][j](q)
6450 * = local_quadrature_points_history[q].old_stress[i][j];
6452 * local_history_strain_values_at_qpoints[i][j](q)
6453 * = local_quadrature_points_history[q].old_strain[i][j];
6455 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
6456 * local_history_stress_values_at_qpoints[i][j]);
6457 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
6458 * history_stress_field[i][j]);
6460 * qpoint_to_dof_matrix.vmult (local_history_strain_fe_values[i][j],
6461 * local_history_strain_values_at_qpoints[i][j]);
6462 * dg_cell->set_dof_values (local_history_strain_fe_values[i][j],
6463 * history_strain_field[i][j]);
6470 * ---------------------------------------------------------------
6474 *
if (refinement_strategy == RefinementStrategy::refine_global)
6479 *
if (cell->is_locally_owned())
6480 * cell->set_refine_flag ();
6484 *
const double refine_fraction_cells = .3,
6485 * coarsen_fraction_cells = .03;
6488 *
const double refine_fraction_cells = .1,
6489 * coarsen_fraction_cells = .3;
6497 * estimated_error_per_cell,
6498 * refine_fraction_cells, coarsen_fraction_cells);
6505 * solution_transfer.prepare_for_coarsening_and_refinement(solution);
6510 *
if (transfer_solution)
6511 * incremental_displacement_transfer.prepare_for_coarsening_and_refinement(incremental_displacement);
6514 * history_stress_field_transfer1(history_dof_handler),
6515 * history_stress_field_transfer2(history_dof_handler);
6516 * history_stress_field_transfer0.prepare_for_coarsening_and_refinement(history_stress_field[0]);
6519 * history_stress_field_transfer1.prepare_for_coarsening_and_refinement(history_stress_field[1]);
6523 * history_stress_field_transfer2.prepare_for_coarsening_and_refinement(history_stress_field[2]);
6527 * history_strain_field_transfer1(history_dof_handler),
6528 * history_strain_field_transfer2(history_dof_handler);
6529 * history_strain_field_transfer0.prepare_for_coarsening_and_refinement(history_strain_field[0]);
6532 * history_strain_field_transfer1.prepare_for_coarsening_and_refinement(history_strain_field[1]);
6536 * history_strain_field_transfer2.prepare_for_coarsening_and_refinement(history_strain_field[2]);
6540 * pcout <<
" Number of active cells: "
6545 * setup_quadrature_point_history ();
6551 * distributed_solution = solution;
6554 * solution_transfer.interpolate(distributed_solution);
6555 * solution = distributed_solution;
6557 *
if (transfer_solution)
6562 * distributed_incremental_displacement = incremental_displacement;
6565 * incremental_displacement_transfer.interpolate(distributed_incremental_displacement);
6566 * incremental_displacement = distributed_incremental_displacement;
6569 * compute_nonlinear_residual(incremental_displacement);
6576 * ---------------------------------------------------
6579 * history_dof_handler.distribute_dofs (history_fe);
6585 * std::vector< std::vector< Vector<double> > >
6586 * distributed_history_stress_field (dim, std::vector<
Vector<double> >(dim));
6587 *
for (
unsigned int i=0; i<dim; ++i)
6588 *
for (
unsigned int j=0; j<dim; ++j)
6590 * distributed_history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6593 * history_stress_field_transfer0.interpolate(history_stress_field[0], distributed_history_stress_field[0]);
6596 * history_stress_field_transfer1.interpolate(history_stress_field[1], distributed_history_stress_field[1]);
6600 * history_stress_field_transfer2.interpolate(history_stress_field[2], distributed_history_stress_field[2]);
6603 * history_stress_field = distributed_history_stress_field;
6610 * std::vector< std::vector< Vector<double> > >
6611 * distributed_history_strain_field (dim, std::vector<
Vector<double> >(dim));
6612 *
for (
unsigned int i=0; i<dim; ++i)
6613 *
for (
unsigned int j=0; j<dim; ++j)
6615 * distributed_history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6618 * history_strain_field_transfer0.interpolate(history_strain_field[0], distributed_history_strain_field[0]);
6621 * history_strain_field_transfer1.interpolate(history_strain_field[1], distributed_history_strain_field[1]);
6625 * history_strain_field_transfer2.interpolate(history_strain_field[2], distributed_history_strain_field[2]);
6628 * history_strain_field = distributed_history_strain_field;
6632 * ---------------------------------------------------------------
6633 * Transfer the history data to the quadrature points of the
new mesh
6634 * In a
final step, we have to get the data back from the now
6635 * interpolated global field to the quadrature points on the
6636 *
new mesh. The following code will
do that:
6643 * history_fe.dofs_per_cell);
6646 * quadrature_formula,
6647 * dof_to_qpoint_matrix);
6648 * cell = dof_handler.begin_active();
6649 * endc = dof_handler.end();
6650 * dg_cell = history_dof_handler.begin_active();
6651 *
for (; cell != endc; ++cell, ++dg_cell)
6652 *
if (cell->is_locally_owned())
6654 * PointHistory<dim> *local_quadrature_points_history
6655 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6656 *
Assert (local_quadrature_points_history >=
6657 * &quadrature_point_history.front(),
6658 * ExcInternalError());
6659 *
Assert (local_quadrature_points_history <
6660 * &quadrature_point_history.back(),
6661 * ExcInternalError());
6662 *
for (
unsigned int i=0; i<dim; ++i)
6663 *
for (
unsigned int j=0; j<dim; ++j)
6665 * dg_cell->get_dof_values (history_stress_field[i][j],
6666 * local_history_stress_fe_values[i][j]);
6667 * dof_to_qpoint_matrix.vmult (local_history_stress_values_at_qpoints[i][j],
6668 * local_history_stress_fe_values[i][j]);
6670 * dg_cell->get_dof_values (history_strain_field[i][j],
6671 * local_history_strain_fe_values[i][j]);
6672 * dof_to_qpoint_matrix.vmult (local_history_strain_values_at_qpoints[i][j],
6673 * local_history_strain_fe_values[i][j]);
6674 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6676 * local_quadrature_points_history[q].old_stress[i][j]
6677 * = local_history_stress_values_at_qpoints[i][j](q);
6679 * local_quadrature_points_history[q].old_strain[i][j]
6680 * = local_history_strain_values_at_qpoints[i][j](q);
6691 * <a name=
"elastoplastic.cc-ElastoPlasticProblemsetup_quadrature_point_history"></a>
6692 * <h4>ElastoPlasticProblem::setup_quadrature_point_history</h4>
6696 * At the beginning of our computations, we needed to
set up
initial values
6697 * of the history variables, such as the existing stresses in the material,
6698 * that we store in each quadrature
point. As mentioned above, we use the
6699 * <code>user_pointer</code>
for this that is available in each cell.
6703 * To put
this into larger perspective, we note that
if we had previously
6704 * available stresses in our model (which we assume
do not exist
for the
6705 * purpose of
this program), then we would need to
interpolate the field of
6706 * preexisting stresses to the quadrature points. Likewise,
if we were to
6707 * simulate elasto-plastic materials with hardening/softening, then we would
6708 * have to store additional history variables like the present yield stress
6709 * of the accumulated plastic strains in each quadrature
6710 * points. Pre-existing hardening or weakening would then be implemented by
6711 * interpolating these variables in the present function as well.
6714 *
template <
int dim>
6715 *
void ElastoPlasticProblem<dim>::setup_quadrature_point_history ()
6719 * What we need to
do here is to
first count how many quadrature points
6720 * are within the responsibility of
this processor. This, of course,
6721 * equals the number of cells that belong to
this processor times the
6722 * number of quadrature points our quadrature formula has on each cell.
6726 * For good measure, we also
set all user pointers of all cells, whether
6727 * ours of not, to the null pointer. This way,
if we ever access the user
6728 * pointer of a cell which we should not have accessed, a segmentation
6729 * fault will let us know that
this should not have happened:
6732 *
unsigned int our_cells = 0;
6736 *
if (cell->is_locally_owned())
6743 * Next, allocate as many quadrature objects as we need. Since the
6744 * <code>resize</code> function does not actually shrink the amount of
6745 * allocated memory
if the requested
new size is smaller than the old
6746 * size, we resort to a trick to
first free all memory, and then
6747 * reallocate it: we declare an empty vector as a temporary variable and
6748 * then
swap the contents of the old vector and
this temporary
6749 * variable. This makes sure that the
6750 * <code>quadrature_point_history</code> is now really empty, and we can
6751 * let the temporary variable that now holds the previous contents of the
6752 * vector go out of scope and be destroyed. In the next step. we can then
6753 * re-allocate as many elements as we need, with the vector
6754 *
default-initializing the <code>PointHistory</code> objects, which
6755 * includes setting the stress variables to zero.
6759 * std::vector<PointHistory<dim> > tmp;
6760 * tmp.swap (quadrature_point_history);
6762 * quadrature_point_history.resize (our_cells *
6763 * quadrature_formula.size());
6767 * Finally
loop over all cells again and
set the user pointers from the
6768 * cells that belong to the present processor to
point to the
first
6769 * quadrature
point objects corresponding to
this cell in the vector of
6773 *
unsigned int history_index = 0;
6777 *
if (cell->is_locally_owned())
6779 * cell->set_user_pointer (&quadrature_point_history[history_index]);
6780 * history_index += quadrature_formula.size();
6785 * At the
end,
for good measure make sure that our count of elements was
6786 * correct and that we have both used up all objects we allocated
6787 * previously, and not
point to any objects beyond the
end of the
6788 * vector. Such defensive programming strategies are
always good checks to
6789 * avoid accidental errors and to guard against future changes to
this
6790 * function that forget to update all uses of a variable at the same
6791 * time. Recall that constructs
using the <code>
Assert</code> macro are
6792 * optimized away in optimized mode, so
do not affect the
run time of
6796 *
Assert (history_index == quadrature_point_history.size(),
6797 * ExcInternalError());
6803 * <a name=
"elastoplastic.cc-ElastoPlasticProblemupdate_quadrature_point_history"></a>
6804 * <h4>ElastoPlasticProblem::update_quadrature_point_history</h4>
6808 * At the
end of each time step, we should have computed an incremental
6809 * displacement update so that the material in its
new configuration
6810 * accommodates
for the difference between the external body and boundary
6811 * forces applied during
this time step minus the forces exerted through
6812 * preexisting
internal stresses. In order to have the preexisting
6813 * stresses available at the next time step, we therefore have to update the
6814 * preexisting stresses with the stresses due to the incremental
6815 * displacement computed during the present time step. Ideally, the
6816 * resulting
sum of
internal stresses would exactly counter all external
6817 * forces. Indeed, a simple experiment can make sure that
this is so:
if we
6818 * choose boundary conditions and body forces to be time
independent, then
6819 * the forcing terms (the sum of external forces and
internal stresses)
6820 * should be exactly zero. If you make
this experiment, you will realize
6821 * from the output of the
norm of the right hand side in each time step that
6822 *
this is almost the
case: it is not exactly zero, since in the
first time
6823 * step the incremental displacement and stress updates were computed
6824 * relative to the undeformed mesh, which was then deformed. In the
second
6825 * time step, we again compute displacement and stress updates, but
this
6826 * time in the deformed mesh -- there, the resulting updates are very small
6827 * but not quite zero. This can be iterated, and in each such iteration the
6828 * residual, i.e. the
norm of the right hand side vector, is reduced;
if one
6829 * makes
this little experiment, one realizes that the
norm of
this residual
6830 * decays exponentially with the number of iterations, and after an
initial
6831 * very rapid decline is reduced by roughly a factor of about 3.5 in each
6832 * iteration (
for one testcase I looked at, other testcases, and other
6833 *
numbers of unknowns change the factor, but not the exponential decay).
6837 * In a sense,
this can then be considered as a quasi-timestepping scheme to
6838 * resolve the nonlinear problem of solving large-deformation elasticity on
6839 * a mesh that is moved along in a Lagrangian manner.
6843 * Another complication is that the existing (old) stresses are defined on
6844 * the old mesh, which we will move around after updating the stresses. If
6845 *
this mesh update involves rotations of the cell, then we need to also
6846 *
rotate the updated stress, since it was computed relative to the
6847 * coordinate system of the old cell.
6851 * Thus, what we need is the following: on each cell which the present
6852 * processor owns, we need to
extract the old stress from the data stored
6853 * with each quadrature
point, compute the stress update, add the two
6854 * together, and then
rotate the result together with the incremental
6855 * rotation computed from the incremental displacement at the present
6856 * quadrature
point. We will detail these steps below:
6859 *
template <
int dim>
6860 *
void ElastoPlasticProblem<dim>::
6861 * update_quadrature_point_history ()
6865 * First,
set up an <code>
FEValues</code>
object by which we will evaluate
6866 * the displacements and the
gradients thereof at the
6867 * quadrature points, together with a vector that will hold
this
6875 *
const unsigned int n_q_points = quadrature_formula.size();
6877 * std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
6883 * Then
loop over all cells and
do the job in the cells that belong to our
6891 * cell = dof_handler.begin_active(),
6892 * endc = dof_handler.end();
6896 *
for (; cell != endc; ++cell)
6897 *
if (cell->is_locally_owned())
6901 * Next, get a pointer to the quadrature
point history data local to
6902 * the present cell, and, as a defensive measure, make sure that
6903 *
this pointer is within the bounds of the global array:
6906 * PointHistory<dim> *local_quadrature_points_history
6907 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6908 *
Assert (local_quadrature_points_history >=
6909 * &quadrature_point_history.front(),
6910 * ExcInternalError());
6911 *
Assert (local_quadrature_points_history <
6912 * &quadrature_point_history.back(),
6913 * ExcInternalError());
6917 * Then initialize the <code>
FEValues</code>
object on the present
6918 * cell, and
extract the strains of the displacement at the
6922 * fe_values.reinit (cell);
6923 * fe_values[displacement].get_function_symmetric_gradients(incremental_displacement,
6924 * incremental_strain_tensor);
6928 * Then
loop over the quadrature points of
this cell:
6931 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6933 * local_quadrature_points_history[q].old_strain +=
6934 * incremental_strain_tensor[q];
6936 * constitutive_law.get_stress_strain_tensor(local_quadrature_points_history[q].old_strain,
6937 * stress_strain_tensor);
6941 * The result of these operations is then written back into
6942 * the original place:
6945 * local_quadrature_points_history[q].old_stress
6946 * = stress_strain_tensor * local_quadrature_points_history[q].old_strain;
6948 * local_quadrature_points_history[q].point
6949 * = fe_values.get_quadrature_points ()[q];
6958 * <a name=
"elastoplastic.cc-PlasticityContactProblemmove_mesh"></a>
6959 * <h4>PlasticityContactProblem::move_mesh</h4>
6963 * The remaining three
functions before we get to <code>
run()</code>
6964 * have to
do with generating output. The following one is an attempt
6965 * at showing the deformed body in its deformed configuration. To
this
6966 *
end,
this function takes a displacement vector field and moves every
6967 * vertex of the (local part) of the mesh by the previously computed
6968 * displacement. We will call
this function with the current
6969 * displacement field before we generate graphical output, and we will
6970 * call it again after generating graphical output with the
negative
6971 * displacement field to undo the changes to the mesh so made.
6975 * The function itself is pretty straightforward. All we have to
do
6976 * is keep track which
vertices we have already touched, as we
6977 * encounter the same
vertices multiple times as we
loop over cells.
6980 *
template <
int dim>
6982 * ElastoPlasticProblem<dim>::
6988 * dof_handler.begin_active();
6989 * cell != dof_handler.end(); ++cell)
6990 *
if (cell->is_locally_owned())
6991 *
for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
6992 *
if (vertex_touched[cell->vertex_index(v)] ==
false)
6994 * vertex_touched[cell->vertex_index(v)] =
true;
6997 *
for (
unsigned int d = 0;
d < dim; ++
d)
6998 * vertex_displacement[d] = displacement(cell->vertex_dof_index(v, d));
7000 * cell->vertex(v) += vertex_displacement;
7009 * <a name=
"elastoplastic.cc-PlasticityContactProblemoutput_results"></a>
7010 * <h4>PlasticityContactProblem::output_results</h4>
7014 * Next is the function we use to actually generate graphical output. The
7015 * function is a bit tedious, but not actually particularly complicated.
7016 * It moves the mesh at the top (and moves it back at the end), then
7017 * computes the contact forces along the contact surface. We can
do
7018 * so (as shown in the accompanying paper) by taking the untreated
7019 * residual vector and identifying which degrees of freedom
7020 * correspond to those with contact by asking whether they have an
7021 * inhomogeneous constraints associated with them. As
always, we need
7022 * to be mindful that we can only write into completely distributed
7023 * vectors (i.e., vectors without ghost elements) but that when we
7024 * want to generate output, we need vectors that
do indeed have
7025 * ghost entries
for all locally relevant degrees of freedom.
7028 *
template <
int dim>
7030 * ElastoPlasticProblem<dim>::output_results (
const std::string &filename_base)
7034 * pcout <<
" Writing graphical output... " << std::flush;
7038 *
const double magnified_factor = 3;
7039 * magnified_solution *= magnified_factor;
7041 * move_mesh(magnified_solution);
7048 *
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
7050 * data_out.add_data_vector(solution,
7051 * std::vector<std::string> (dim,
"displacement"),
7055 * std::vector<std::string> solution_names;
7060 * solution_names.push_back (
"displacement");
7063 * solution_names.push_back (
"x_displacement");
7064 * solution_names.push_back (
"y_displacement");
7067 * solution_names.push_back (
"x_displacement");
7068 * solution_names.push_back (
"y_displacement");
7069 * solution_names.push_back (
"z_displacement");
7075 * data_out.add_data_vector (solution, solution_names);
7080 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
7082 * data_out.add_data_vector(subdomain,
"subdomain");
7085 * data_out.add_data_vector(fraction_of_plastic_q_points_per_cell,
7086 *
"fraction_of_plastic_q_points");
7089 * data_out.build_patches();
7093 * In the remainder of the function, we generate one VTU file on
7094 * every processor, indexed by the subdomain
id of
this processor.
7095 * On the
first processor, we then also create a <code>.pvtu</code>
7096 * file that indexes <i>all</i> of the VTU files so that the entire
7097 *
set of output files can be read at once. These <code>.pvtu</code>
7098 * are used by Paraview to describe an entire
parallel computation
's
7099 * output files. We then do the same again for the competitor of
7100 * Paraview, the Visit visualization program, by creating a matching
7101 * <code>.visit</code> file.
7104 * const std::string filename =
7105 * (output_dir + filename_base + "-"
7106 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7108 * std::ofstream output_vtu((filename + ".vtu").c_str());
7109 * data_out.write_vtu(output_vtu);
7110 * pcout << output_dir + filename_base << ".pvtu" << std::endl;
7113 * if (this_mpi_process == 0)
7115 * std::vector<std::string> filenames;
7116 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7117 * filenames.push_back(filename_base + "-" +
7118 * Utilities::int_to_string(i, 4) +
7121 * std::ofstream pvtu_master_output((output_dir + filename_base + ".pvtu").c_str());
7122 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7124 * std::ofstream visit_master_output((output_dir + filename_base + ".visit").c_str());
7125 * data_out.write_pvtu_record(visit_master_output, filenames);
7129 * produce eps files for mesh illustration
7132 * std::ofstream output_eps((filename + ".eps").c_str());
7134 * grid_out.write_eps(triangulation, output_eps);
7139 * Extrapolate the stresses from Gauss point to the nodes
7142 * SymmetricTensor<2, dim> stress_at_qpoint;
7144 * FE_DGQ<dim> history_fe (1);
7145 * DoFHandler<dim> history_dof_handler (triangulation);
7146 * history_dof_handler.distribute_dofs (history_fe);
7147 * std::vector< std::vector< Vector<double> > >
7148 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7149 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7150 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7151 * for (unsigned int i=0; i<dim; ++i)
7152 * for (unsigned int j=0; j<dim; ++j)
7154 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7155 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7156 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7159 * Vector<double> VM_stress_field (history_dof_handler.n_dofs()),
7160 * local_VM_stress_values_at_qpoints (quadrature_formula.size()),
7161 * local_VM_stress_fe_values (history_fe.dofs_per_cell);
7163 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7164 * quadrature_formula.size());
7165 * FETools::compute_projection_from_quadrature_points_matrix
7167 * quadrature_formula, quadrature_formula,
7168 * qpoint_to_dof_matrix);
7170 * typename DoFHandler<dim>::active_cell_iterator
7171 * cell = dof_handler.begin_active(),
7172 * endc = dof_handler.end(),
7173 * dg_cell = history_dof_handler.begin_active();
7175 * const FEValuesExtractors::Vector displacement(0);
7177 * for (; cell!=endc; ++cell, ++dg_cell)
7178 * if (cell->is_locally_owned())
7180 * PointHistory<dim> *local_quadrature_points_history
7181 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7182 * Assert (local_quadrature_points_history >=
7183 * &quadrature_point_history.front(),
7184 * ExcInternalError());
7185 * Assert (local_quadrature_points_history <
7186 * &quadrature_point_history.back(),
7187 * ExcInternalError());
7191 * Then loop over the quadrature points of this cell:
7194 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7196 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7198 * for (unsigned int i=0; i<dim; ++i)
7199 * for (unsigned int j=i; j<dim; ++j)
7201 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7204 * local_VM_stress_values_at_qpoints(q) = Evaluation::get_von_Mises_stress(stress_at_qpoint);
7209 * for (unsigned int i=0; i<dim; ++i)
7210 * for (unsigned int j=i; j<dim; ++j)
7212 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7213 * local_history_stress_values_at_qpoints[i][j]);
7214 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7215 * history_stress_field[i][j]);
7218 * qpoint_to_dof_matrix.vmult (local_VM_stress_fe_values,
7219 * local_VM_stress_values_at_qpoints);
7220 * dg_cell->set_dof_values (local_VM_stress_fe_values,
7228 * Save stresses on nodes by nodal averaging
7229 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7230 * in order to compute stresses on nodes (by applying nodal averaging)
7231 * Therefore, each vertex has one degree of freedom
7234 * FE_Q<dim> fe_1 (1);
7235 * DoFHandler<dim> dof_handler_1 (triangulation);
7236 * dof_handler_1.distribute_dofs (fe_1);
7238 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7239 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7241 * std::vector< std::vector< Vector<double> > >
7242 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7243 * for (unsigned int i=0; i<dim; ++i)
7244 * for (unsigned int j=0; j<dim; ++j)
7246 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7249 * Vector<double> VM_stress_on_vertices (dof_handler_1.n_dofs()),
7250 * counter_on_vertices (dof_handler_1.n_dofs());
7251 * VM_stress_on_vertices = 0;
7252 * counter_on_vertices = 0;
7254 * cell = dof_handler.begin_active();
7255 * dg_cell = history_dof_handler.begin_active();
7256 * typename DoFHandler<dim>::active_cell_iterator
7257 * cell_1 = dof_handler_1.begin_active();
7258 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7259 * if (cell->is_locally_owned())
7261 * dg_cell->get_dof_values (VM_stress_field,
7262 * local_VM_stress_fe_values);
7264 * for (unsigned int i=0; i<dim; ++i)
7265 * for (unsigned int j=0; j<dim; ++j)
7267 * dg_cell->get_dof_values (history_stress_field[i][j],
7268 * local_history_stress_fe_values[i][j]);
7271 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7273 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7278 * Point<dim> point1, point2;
7279 * point1 = cell_1->vertex(v);
7280 * point2 = dg_cell->vertex(v);
7281 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7288 * counter_on_vertices (dof_1_vertex) += 1;
7290 * VM_stress_on_vertices (dof_1_vertex) += local_VM_stress_fe_values (v);
7292 * for (unsigned int i=0; i<dim; ++i)
7293 * for (unsigned int j=0; j<dim; ++j)
7295 * history_stress_on_vertices[i][j](dof_1_vertex) +=
7296 * local_history_stress_fe_values[i][j](v);
7302 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7304 * VM_stress_on_vertices(id) /= counter_on_vertices(id);
7306 * for (unsigned int i=0; i<dim; ++i)
7307 * for (unsigned int j=0; j<dim; ++j)
7309 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7315 * Save figures of stresses
7318 * if (show_stresses)
7321 * DataOut<dim> data_out;
7322 * data_out.attach_dof_handler (history_dof_handler);
7325 * data_out.add_data_vector (history_stress_field[0][0], "stress_xx");
7326 * data_out.add_data_vector (history_stress_field[1][1], "stress_yy");
7327 * data_out.add_data_vector (history_stress_field[0][1], "stress_xy");
7328 * data_out.add_data_vector (VM_stress_field, "Von_Mises_stress");
7332 * data_out.add_data_vector (history_stress_field[0][2], "stress_xz");
7333 * data_out.add_data_vector (history_stress_field[1][2], "stress_yz");
7334 * data_out.add_data_vector (history_stress_field[2][2], "stress_zz");
7337 * data_out.build_patches ();
7339 * const std::string filename_base_stress = ("stress-" + filename_base);
7341 * const std::string filename =
7342 * (output_dir + filename_base_stress + "-"
7343 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7345 * std::ofstream output_vtu((filename + ".vtu").c_str());
7346 * data_out.write_vtu(output_vtu);
7347 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7349 * if (this_mpi_process == 0)
7351 * std::vector<std::string> filenames;
7352 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7353 * filenames.push_back(filename_base_stress + "-" +
7354 * Utilities::int_to_string(i, 4) +
7357 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7358 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7360 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7361 * data_out.write_pvtu_record(visit_master_output, filenames);
7368 * DataOut<dim> data_out;
7369 * data_out.attach_dof_handler (dof_handler_1);
7372 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_xx_averaged");
7373 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_yy_averaged");
7374 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_xy_averaged");
7375 * data_out.add_data_vector (VM_stress_on_vertices, "Von_Mises_stress_averaged");
7379 * data_out.add_data_vector (history_stress_on_vertices[0][2], "stress_xz_averaged");
7380 * data_out.add_data_vector (history_stress_on_vertices[1][2], "stress_yz_averaged");
7381 * data_out.add_data_vector (history_stress_on_vertices[2][2], "stress_zz_averaged");
7384 * data_out.build_patches ();
7386 * const std::string filename_base_stress = ("averaged-stress-" + filename_base);
7388 * const std::string filename =
7389 * (output_dir + filename_base_stress + "-"
7390 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7392 * std::ofstream output_vtu((filename + ".vtu").c_str());
7393 * data_out.write_vtu(output_vtu);
7394 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7396 * if (this_mpi_process == 0)
7398 * std::vector<std::string> filenames;
7399 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7400 * filenames.push_back(filename_base_stress + "-" +
7401 * Utilities::int_to_string(i, 4) +
7404 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7405 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7407 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7408 * data_out.write_pvtu_record(visit_master_output, filenames);
7415 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7423 * magnified_solution *= -1;
7424 * move_mesh(magnified_solution);
7431 * if (base_mesh == "Timoshenko beam")
7433 * const double length = .48,
7436 * Point<dim> intersted_point(length, -depth/2);
7437 * Point<dim> vertex_displacement;
7438 * bool vertex_found = false;
7440 * for (typename DoFHandler<dim>::active_cell_iterator cell =
7441 * dof_handler.begin_active();
7442 * cell != dof_handler.end(); ++cell)
7443 * if (cell->is_locally_owned() && !vertex_found)
7444 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7445 * if ( std::fabs(cell->vertex(v)[0] - intersted_point[0])<1e-6 &&
7446 * std::fabs(cell->vertex(v)[1] - intersted_point[1])<1e-6)
7448 * vertex_found = true;
7450 * for (unsigned int d = 0; d < dim; ++d)
7451 * vertex_displacement[d] = solution(cell->vertex_dof_index(v, d));
7456 * pcout << " Number of active cells: "
7457 * << triangulation.n_global_active_cells() << std::endl
7458 * << " Number of degrees of freedom: " << dof_handler.n_dofs()
7461 * AssertThrow(vertex_found, ExcInternalError());
7462 * std::cout << "Displacement at the point (" << intersted_point[0]
7463 * << ", " << intersted_point[1] << ") is "
7464 * << "(" << vertex_displacement[0]
7465 * << ", " << vertex_displacement[1] << ").\n";
7467 * Vector<double> vertex_exact_displacement(dim);
7468 * EquationData::IncrementalBoundaryValues<dim> incremental_boundary_values(present_time, end_time);
7469 * incremental_boundary_values.vector_value (intersted_point, vertex_exact_displacement);
7471 * std::cout << "Exact displacement at the point (" << intersted_point[0]
7472 * << ", " << intersted_point[1] << ") is "
7473 * << "(" << vertex_exact_displacement[0]
7474 * << ", " << vertex_exact_displacement[1] << ").\n\n";
7477 * else if (base_mesh == "Thick_tube_internal_pressure")
7479 * const double pressure (0.6*2.4e8),
7480 * inner_radius (.1);
7483 * const double pressure (1.94e8),
7484 * inner_radius (.1);
7492 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7493 * 3d and plane strain
7496 * const double mu (e_modulus / (2 * (1 + nu)));
7498 * const Point<dim> point_A(inner_radius, 0.);
7499 * Vector<double> disp_A(dim);
7503 * make a non-parallel copy of solution
7506 * Vector<double> copy_solution(solution);
7508 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
7510 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
7512 * table_results.add_value("time step", timestep_no);
7513 * table_results.add_value("Cells", triangulation.n_global_active_cells());
7514 * table_results.add_value("DoFs", dof_handler.n_dofs());
7515 * table_results.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7516 * table_results.add_value("4*mu*u_A/(sigma_0*a)", 4*mu*disp_A(0)/(sigma_0*inner_radius));
7520 * Compute stresses in the POLAR coordinates, 1- save it on Gauss points,
7521 * 2- extrapolate them to nodes and taking their avarages (nodal avaraging)
7524 * AssertThrow (dim == 2, ExcNotImplemented());
7528 * we define a rotation matrix to be able to transform the stress
7529 * from the Cartesian coordinate to the polar coordinate
7532 * Tensor<2, dim> rotation_matrix; // [cos sin; -sin cos] , sigma_r = rot * sigma * rot^T
7534 * FEValues<dim> fe_values (fe, quadrature_formula, update_quadrature_points |
7535 * update_values | update_gradients);
7537 * const unsigned int n_q_points = quadrature_formula.size();
7539 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
7540 * SymmetricTensor<4, dim> stress_strain_tensor;
7541 * Tensor<2, dim> stress_at_qpoint;
7543 * FE_DGQ<dim> history_fe (1);
7544 * DoFHandler<dim> history_dof_handler (triangulation);
7545 * history_dof_handler.distribute_dofs (history_fe);
7546 * std::vector< std::vector< Vector<double> > >
7547 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7548 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7549 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7550 * for (unsigned int i=0; i<dim; ++i)
7551 * for (unsigned int j=0; j<dim; ++j)
7553 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7554 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7555 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7558 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7559 * quadrature_formula.size());
7560 * FETools::compute_projection_from_quadrature_points_matrix
7562 * quadrature_formula, quadrature_formula,
7563 * qpoint_to_dof_matrix);
7565 * typename DoFHandler<dim>::active_cell_iterator
7566 * cell = dof_handler.begin_active(),
7567 * endc = dof_handler.end(),
7568 * dg_cell = history_dof_handler.begin_active();
7570 * const FEValuesExtractors::Vector displacement(0);
7572 * for (; cell!=endc; ++cell, ++dg_cell)
7573 * if (cell->is_locally_owned())
7575 * PointHistory<dim> *local_quadrature_points_history
7576 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7577 * Assert (local_quadrature_points_history >=
7578 * &quadrature_point_history.front(),
7579 * ExcInternalError());
7580 * Assert (local_quadrature_points_history <
7581 * &quadrature_point_history.back(),
7582 * ExcInternalError());
7586 * Then loop over the quadrature points of this cell:
7589 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7591 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7595 * transform the stress from the Cartesian coordinate to the polar coordinate
7598 * const Point<dim> point = local_quadrature_points_history[q].point;
7599 * const double theta = std::atan2(point(1),point(0));
7606 * rotation_matrix[0][0] = std::cos(theta);
7607 * rotation_matrix[0][1] = std::sin(theta);
7608 * rotation_matrix[1][0] = -std::sin(theta);
7609 * rotation_matrix[1][1] = std::cos(theta);
7613 * stress in polar coordinate
7616 * stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix);
7618 * for (unsigned int i=0; i<dim; ++i)
7619 * for (unsigned int j=i; j<dim; ++j)
7621 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7627 * for (unsigned int i=0; i<dim; ++i)
7628 * for (unsigned int j=i; j<dim; ++j)
7630 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7631 * local_history_stress_values_at_qpoints[i][j]);
7632 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7633 * history_stress_field[i][j]);
7639 * DataOut<dim> data_out;
7640 * data_out.attach_dof_handler (history_dof_handler);
7643 * data_out.add_data_vector (history_stress_field[0][0], "stress_rr");
7644 * data_out.add_data_vector (history_stress_field[1][1], "stress_tt");
7645 * data_out.add_data_vector (history_stress_field[0][1], "stress_rt");
7647 * data_out.build_patches ();
7649 * const std::string filename_base_stress = ("stress-polar-" + filename_base);
7651 * const std::string filename =
7652 * (output_dir + filename_base_stress + "-"
7653 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7655 * std::ofstream output_vtu((filename + ".vtu").c_str());
7656 * data_out.write_vtu(output_vtu);
7657 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7659 * if (this_mpi_process == 0)
7661 * std::vector<std::string> filenames;
7662 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7663 * filenames.push_back(filename_base_stress + "-" +
7664 * Utilities::int_to_string(i, 4) +
7667 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7668 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7670 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7671 * data_out.write_pvtu_record(visit_master_output, filenames);
7679 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7680 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7681 * in order to compute stresses on nodes (by applying nodal averaging)
7682 * Therefore, each vertex has one degree of freedom
7685 * FE_Q<dim> fe_1 (1);
7686 * DoFHandler<dim> dof_handler_1 (triangulation);
7687 * dof_handler_1.distribute_dofs (fe_1);
7689 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7690 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7692 * std::vector< std::vector< Vector<double> > >
7693 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7694 * for (unsigned int i=0; i<dim; ++i)
7695 * for (unsigned int j=0; j<dim; ++j)
7697 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7700 * Vector<double> counter_on_vertices (dof_handler_1.n_dofs());
7701 * counter_on_vertices = 0;
7703 * cell = dof_handler.begin_active();
7704 * dg_cell = history_dof_handler.begin_active();
7705 * typename DoFHandler<dim>::active_cell_iterator
7706 * cell_1 = dof_handler_1.begin_active();
7707 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7708 * if (cell->is_locally_owned())
7711 * for (unsigned int i=0; i<dim; ++i)
7712 * for (unsigned int j=0; j<dim; ++j)
7714 * dg_cell->get_dof_values (history_stress_field[i][j],
7715 * local_history_stress_fe_values[i][j]);
7718 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7720 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7725 * Point<dim> point1, point2;
7726 * point1 = cell_1->vertex(v);
7727 * point2 = dg_cell->vertex(v);
7728 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7735 * counter_on_vertices (dof_1_vertex) += 1;
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](dof_1_vertex) +=
7741 * local_history_stress_fe_values[i][j](v);
7747 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7749 * for (unsigned int i=0; i<dim; ++i)
7750 * for (unsigned int j=0; j<dim; ++j)
7752 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7758 * DataOut<dim> data_out;
7759 * data_out.attach_dof_handler (dof_handler_1);
7762 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_rr_averaged");
7763 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_tt_averaged");
7764 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_rt_averaged");
7766 * data_out.build_patches ();
7768 * const std::string filename_base_stress = ("averaged-stress-polar-" + filename_base);
7770 * const std::string filename =
7771 * (output_dir + filename_base_stress + "-"
7772 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7774 * std::ofstream output_vtu((filename + ".vtu").c_str());
7775 * data_out.write_vtu(output_vtu);
7776 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7778 * if (this_mpi_process == 0)
7780 * std::vector<std::string> filenames;
7781 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7782 * filenames.push_back(filename_base_stress + "-" +
7783 * Utilities::int_to_string(i, 4) +
7786 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7787 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7789 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7790 * data_out.write_pvtu_record(visit_master_output, filenames);
7797 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7803 * if ( std::abs( (present_time/end_time)*(pressure/sigma_0) - 0.6 ) <
7804 * .501*(present_timestep/end_time)*(pressure/sigma_0) )
7809 * table_results_2: presenting the stress_rr and stress_tt on the nodes of bottom edge
7812 * const unsigned int face_id = 3;
7814 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
7816 * bool evaluation_face_found = false;
7818 * typename DoFHandler<dim>::active_cell_iterator
7819 * cell = dof_handler.begin_active(),
7820 * endc = dof_handler.end(),
7821 * cell_1 = dof_handler_1.begin_active();
7822 * for (; cell!=endc; ++cell, ++cell_1)
7823 * if (cell->is_locally_owned())
7825 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7827 * if (cell->face(face)->at_boundary()
7829 * cell->face(face)->boundary_id() == face_id)
7831 * if (!evaluation_face_found)
7833 * evaluation_face_found = true;
7837 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7839 * types::global_dof_index dof_1_vertex =
7840 * cell_1->face(face)->vertex_dof_index(v, 0);
7841 * if (!vertices_found[dof_1_vertex])
7844 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
7846 * table_results_2.add_value("x coordinate", vertex_coordinate[0]);
7847 * table_results_2.add_value("stress_rr", history_stress_on_vertices[0][0](dof_1_vertex));
7848 * table_results_2.add_value("stress_tt", history_stress_on_vertices[1][1](dof_1_vertex));
7849 * table_results_2.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7851 * vertices_found[dof_1_vertex] = true;
7860 * AssertThrow(evaluation_face_found, ExcInternalError());
7864 * table_results_3: presenting the mean stress_rr of the nodes on the inner radius
7867 * const unsigned int face_id_2 = 0;
7869 * Tensor<2, dim> stress_node,
7870 * mean_stress_polar;
7871 * mean_stress_polar = 0;
7873 * std::vector<bool> vertices_found_2 (dof_handler_1.n_dofs(), false);
7874 * unsigned int no_vertices_found = 0;
7876 * evaluation_face_found = false;
7878 * cell = dof_handler.begin_active(),
7879 * endc = dof_handler.end(),
7880 * cell_1 = dof_handler_1.begin_active();
7881 * for (; cell!=endc; ++cell, ++cell_1)
7882 * if (cell->is_locally_owned())
7884 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7886 * if (cell->face(face)->at_boundary()
7888 * cell->face(face)->boundary_id() == face_id_2)
7890 * if (!evaluation_face_found)
7892 * evaluation_face_found = true;
7896 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7898 * types::global_dof_index dof_1_vertex =
7899 * cell_1->face(face)->vertex_dof_index(v, 0);
7900 * if (!vertices_found_2[dof_1_vertex])
7902 * for (unsigned int ir=0; ir<dim; ++ir)
7903 * for (unsigned int ic=0; ic<dim; ++ic)
7904 * stress_node[ir][ic] = history_stress_on_vertices[ir][ic](dof_1_vertex);
7906 * mean_stress_polar += stress_node;
7908 * vertices_found_2[dof_1_vertex] = true;
7909 * ++no_vertices_found;
7918 * AssertThrow(evaluation_face_found, ExcInternalError());
7920 * mean_stress_polar /= no_vertices_found;
7922 * table_results_3.add_value("time step", timestep_no);
7923 * table_results_3.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7924 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
7925 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
7926 * table_results_3.add_value("radius", inner_radius);
7927 * table_results_3.add_value("mean stress_rr", mean_stress_polar[0][0]);
7928 * table_results_3.add_value("mean stress_tt", mean_stress_polar[1][1]);
7935 * else if (base_mesh == "Perforated_strip_tension")
7937 * const double imposed_displacement (0.00055),
7938 * inner_radius (0.05);
7943 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7944 * 3d and plane strain
7948 * table_results: Demonstrates the result of displacement at the top left corner versus imposed tension
7953 * const Point<dim> point_C(0., height);
7954 * Vector<double> disp_C(dim);
7958 * make a non-parallel copy of solution
7961 * Vector<double> copy_solution(solution);
7963 * typename Evaluation::PointValuesEvaluation<dim>::
7964 * PointValuesEvaluation point_values_evaluation(point_C);
7966 * point_values_evaluation.compute (dof_handler, copy_solution, disp_C);
7968 * table_results.add_value("time step", timestep_no);
7969 * table_results.add_value("Cells", triangulation.n_global_active_cells());
7970 * table_results.add_value("DoFs", dof_handler.n_dofs());
7971 * table_results.add_value("4*mu*u_C/(sigma_0*r)", 4*mu*disp_C(1)/(sigma_0*inner_radius));
7977 * compute average sigma_yy on the bottom edge
7980 * double stress_yy_av;
7983 * const unsigned int face_id = 1;
7985 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
7986 * unsigned int no_vertices_in_face = 0;
7988 * bool evaluation_face_found = false;
7990 * typename DoFHandler<dim>::active_cell_iterator
7991 * cell = dof_handler.begin_active(),
7992 * endc = dof_handler.end(),
7993 * cell_1 = dof_handler_1.begin_active();
7994 * for (; cell!=endc; ++cell, ++cell_1)
7995 * if (cell->is_locally_owned())
7997 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7999 * if (cell->face(face)->at_boundary()
8001 * cell->face(face)->boundary_id() == face_id)
8003 * if (!evaluation_face_found)
8005 * evaluation_face_found = true;
8009 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
8011 * types::global_dof_index dof_1_vertex =
8012 * cell_1->face(face)->vertex_dof_index(v, 0);
8013 * if (!vertices_found[dof_1_vertex])
8015 * stress_yy_av += history_stress_on_vertices[1][1](dof_1_vertex);
8016 * ++no_vertices_in_face;
8018 * vertices_found[dof_1_vertex] = true;
8027 * AssertThrow(evaluation_face_found, ExcInternalError());
8029 * stress_yy_av /= no_vertices_in_face;
8035 * table_results_2: Demonstrate the stress_yy on the nodes of bottom edge
8039 * if ( std::abs( (stress_yy_av/sigma_0) - .91 ) < .2 )
8042 * if ( (timestep_no) % 19 == 0 )
8049 * const unsigned int face_id = 1;
8051 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
8053 * bool evaluation_face_found = false;
8055 * typename DoFHandler<dim>::active_cell_iterator
8056 * cell = dof_handler.begin_active(),
8057 * endc = dof_handler.end(),
8058 * cell_1 = dof_handler_1.begin_active();
8059 * for (; cell!=endc; ++cell, ++cell_1)
8060 * if (cell->is_locally_owned())
8062 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
8064 * if (cell->face(face)->at_boundary()
8066 * cell->face(face)->boundary_id() == face_id)
8068 * if (!evaluation_face_found)
8070 * evaluation_face_found = true;
8074 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
8076 * types::global_dof_index dof_1_vertex =
8077 * cell_1->face(face)->vertex_dof_index(v, 0);
8079 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
8081 * if (!vertices_found[dof_1_vertex] && std::abs(vertex_coordinate[2])<1.e-8)
8083 * table_results_2.add_value("x", vertex_coordinate[0]);
8084 * table_results_2.add_value("x/r", vertex_coordinate[0]/inner_radius);
8085 * table_results_2.add_value("stress_xx/sigma_0", history_stress_on_vertices[0][0](dof_1_vertex)/sigma_0);
8086 * table_results_2.add_value("stress_yy/sigma_0", history_stress_on_vertices[1][1](dof_1_vertex)/sigma_0);
8087 * table_results_2.add_value("stress_yy_av/sigma_0", stress_yy_av/sigma_0);
8088 * table_results_2.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8090 * vertices_found[dof_1_vertex] = true;
8099 * AssertThrow(evaluation_face_found, ExcInternalError());
8105 * table_results_3: Demonstrate the Stress_mean (average tensile stress)
8106 * on the bottom edge versus epsilon_yy on the bottom left corner
8110 * double strain_yy_A = 0.;
8114 * compute strain_yy_A
8115 * Since the point A is the node on the bottom left corner,
8116 * we need to work just with one element
8120 * const Point<dim> point_A(inner_radius, 0, 0);
8122 * Vector<double> local_strain_yy_values_at_qpoints (quadrature_formula.size()),
8123 * local_strain_yy_fe_values (history_fe.dofs_per_cell);
8125 * SymmetricTensor<2, dim> strain_at_qpoint;
8127 * typename DoFHandler<dim>::active_cell_iterator
8128 * cell = dof_handler.begin_active(),
8129 * endc = dof_handler.end(),
8130 * dg_cell = history_dof_handler.begin_active();
8132 * bool cell_found = false;
8134 * for (; cell!=endc; ++cell, ++dg_cell)
8135 * if (cell->is_locally_owned() && !cell_found)
8137 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
8138 * if ( std::fabs(cell->vertex(v)[0] - point_A[0])<1e-6 &&
8139 * std::fabs(cell->vertex(v)[1] - point_A[1])<1e-6 &&
8140 * std::fabs(cell->vertex(v)[2] - point_A[2])<1e-6)
8142 * PointHistory<dim> *local_quadrature_points_history
8143 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8144 * Assert (local_quadrature_points_history >=
8145 * &quadrature_point_history.front(),
8146 * ExcInternalError());
8147 * Assert (local_quadrature_points_history <
8148 * &quadrature_point_history.back(),
8149 * ExcInternalError());
8153 * Then loop over the quadrature points of this cell:
8156 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8158 * strain_at_qpoint = local_quadrature_points_history[q].old_strain;
8160 * local_strain_yy_values_at_qpoints(q) = strain_at_qpoint[1][1];
8163 * qpoint_to_dof_matrix.vmult (local_strain_yy_fe_values,
8164 * local_strain_yy_values_at_qpoints);
8166 * strain_yy_A = local_strain_yy_fe_values (v);
8168 * cell_found = true;
8176 * table_results_3.add_value("time step", timestep_no);
8177 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
8178 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
8179 * table_results_3.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8180 * table_results_3.add_value("mean_tensile_stress/sigma_0", stress_yy_av/sigma_0);
8181 * table_results_3.add_value("E*strain_yy-A/sigma_0", e_modulus*strain_yy_A/sigma_0);
8186 * if (std::abs(present_time-end_time) < 1.e-7)
8188 * table_results_2.set_precision("Imposed u_y", 6);
8189 * table_results_3.set_precision("Imposed u_y", 6);
8193 * else if (base_mesh == "Cantiliver_beam_3d")
8195 * const double pressure (6e6),
8201 * table_results: Demonstrates the result of displacement at the top front point, Point A
8207 * Quantity of interest:
8208 * displacement at Point A (x=0, y=height/2, z=length)
8214 * const Point<dim> point_A(0, height/2, length);
8215 * Vector<double> disp_A(dim);
8219 * make a non-parallel copy of solution
8222 * Vector<double> copy_solution(solution);
8224 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
8226 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
8228 * table_results.add_value("time step", timestep_no);
8229 * table_results.add_value("Cells", triangulation.n_global_active_cells());
8230 * table_results.add_value("DoFs", dof_handler.n_dofs());
8231 * table_results.add_value("pressure", pressure*present_time/end_time);
8232 * table_results.add_value("u_A", disp_A(1));
8238 * demonstrate the location and maximum von-Mises stress in the
8239 * specified domain close to the clamped face, z = 0
8240 * top domain: height/2 - thickness_flange <= y <= height/2
8241 * 0 <= z <= 2 * thickness_flange
8242 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
8243 * 0 <= z <= 2 * thickness_flange
8249 * double VM_stress_max (0);
8250 * Point<dim> point_max;
8252 * SymmetricTensor<2, dim> stress_at_qpoint;
8254 * typename DoFHandler<dim>::active_cell_iterator
8255 * cell = dof_handler.begin_active(),
8256 * endc = dof_handler.end();
8258 * const FEValuesExtractors::Vector displacement(0);
8260 * for (; cell!=endc; ++cell)
8261 * if (cell->is_locally_owned())
8263 * PointHistory<dim> *local_quadrature_points_history
8264 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8265 * Assert (local_quadrature_points_history >=
8266 * &quadrature_point_history.front(),
8267 * ExcInternalError());
8268 * Assert (local_quadrature_points_history <
8269 * &quadrature_point_history.back(),
8270 * ExcInternalError());
8274 * Then loop over the quadrature points of this cell:
8277 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8279 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
8281 * const double VM_stress = Evaluation::get_von_Mises_stress(stress_at_qpoint);
8282 * if (VM_stress > VM_stress_max)
8284 * VM_stress_max = VM_stress;
8285 * point_max = local_quadrature_points_history[q].point;
8291 * table_results.add_value("maximum von_Mises stress", VM_stress_max);
8292 * table_results.add_value("x", point_max[0]);
8293 * table_results.add_value("y", point_max[1]);
8294 * table_results.add_value("z", point_max[2]);
8307 * <a name="elastoplastic.cc-PlasticityContactProblemrun"></a>
8308 * <h4>PlasticityContactProblem::run</h4>
8312 * As in all other tutorial programs, the <code>run()</code> function contains
8313 * the overall logic. There is not very much to it here: in essence, it
8314 * performs the loops over all mesh refinement cycles, and within each, hands
8315 * things over to the Newton solver in <code>solve_newton()</code> on the
8316 * current mesh and calls the function that creates graphical output for
8317 * the so-computed solution. It then outputs some statistics concerning both
8318 * run times and memory consumption that has been collected over the course of
8319 * computations on this mesh.
8322 * template <int dim>
8324 * ElastoPlasticProblem<dim>::run ()
8326 * computing_timer.reset();
8329 * present_timestep = 1;
8337 * ----------------------------------------------------------------
8338 * base_mesh == "Thick_tube_internal_pressure"
8342 * const Point<dim> center(0, 0);
8343 * const double inner_radius = .1,
8344 * outer_radius = .2;
8346 * const SphericalManifold<dim> inner_boundary_description(center, inner_radius);
8347 * triangulation.set_manifold (0, inner_boundary_description);
8349 * const SphericalManifold<dim> outer_boundary_description(center, outer_radius);
8350 * triangulation.set_manifold (1, outer_boundary_description);
8354 * ----------------------------------------------------------------
8355 * base_mesh == "Perforated_strip_tension"
8359 * const double inner_radius = 0.05;
8361 * const CylinderBoundary<dim> inner_boundary_description(inner_radius, 2);
8362 * triangulation.set_manifold (10, inner_boundary_description);
8366 * ----------------------------------------------------------------
8372 * setup_quadrature_point_history ();
8374 * while (present_time < end_time)
8376 * present_time += present_timestep;
8379 * if (present_time > end_time)
8381 * present_timestep -= (present_time - end_time);
8382 * present_time = end_time;
8384 * pcout << std::endl;
8385 * pcout << "Time step " << timestep_no << " at time " << present_time
8388 * relative_error = max_relative_error * 10;
8389 * current_refinement_cycle = 0;
8396 * ------------------------ Refinement based on the relative error -------------------------------
8402 * while (relative_error >= max_relative_error)
8407 * if ( (timestep_no > 1) && (current_refinement_cycle>0) && (relative_error >= max_relative_error) )
8409 * pcout << "The relative error, " << relative_error
8410 * << " , is still more than maximum relative error, "
8411 * << max_relative_error << ", but we move to the next increment.\n";
8412 * relative_error = .1 * max_relative_error;
8415 * if (relative_error >= max_relative_error)
8417 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8418 * ++current_refinement_cycle;
8426 * ------------------------ Refinement based on the number of refinement --------------------------
8430 * bool continue_loop = true;
8431 * while (continue_loop)
8436 * if ( (timestep_no == 1) && (current_refinement_cycle < 1) )
8438 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8439 * ++current_refinement_cycle;
8443 * continue_loop = false;
8451 * -------------------------------------------------------------------------------------------------
8457 * solution += incremental_displacement;
8459 * update_quadrature_point_history ();
8461 * output_results((std::string("solution-") +
8462 * Utilities::int_to_string(timestep_no, 4)).c_str());
8464 * computing_timer.print_summary();
8465 * computing_timer.reset();
8467 * Utilities::System::MemoryStats stats;
8468 * Utilities::System::get_memory_stats(stats);
8469 * pcout << "Peak virtual memory used, resident in kB: " << stats.VmSize << " "
8470 * << stats.VmRSS << std::endl;
8473 * if (std::abs(present_time-end_time) < 1.e-7)
8475 * const std::string filename = (output_dir + "Results");
8477 * std::ofstream output_txt((filename + ".txt").c_str());
8479 * pcout << std::endl;
8480 * table_results.write_text(output_txt);
8481 * pcout << std::endl;
8482 * table_results_2.write_text(output_txt);
8483 * pcout << std::endl;
8484 * table_results_3.write_text(output_txt);
8485 * pcout << std::endl;
8490 * if (base_mesh == "Thick_tube_internal_pressure")
8492 * triangulation.reset_manifold (0);
8493 * triangulation.reset_manifold (1);
8495 * else if (base_mesh == "Perforated_strip_tension")
8497 * triangulation.reset_manifold (10);
8506 * <a name="elastoplastic.cc-Thecodemaincodefunction"></a>
8507 * <h3>The <code>main</code> function</h3>
8511 * There really isn't much to the <code>main()</code> function. It looks
8515 *
int main (
int argc,
char *argv[])
8517 *
using namespace dealii;
8518 *
using namespace ElastoPlastic;
8524 *
const int dim = 3;
8525 * ElastoPlasticProblem<dim>::declare_parameters(prm);
8528 * std::cerr <<
"*** Call this program as <./elastoplastic input.prm>" << std::endl;
8532 * prm.parse_input(argv[1]);
8535 * ElastoPlasticProblem<dim> problem(prm);
8539 *
catch (std::exception &exc)
8541 * std::cerr << std::endl << std::endl
8542 * <<
"----------------------------------------------------"
8544 * std::cerr <<
"Exception on processing: " << std::endl
8545 * << exc.what() << std::endl
8546 * <<
"Aborting!" << std::endl
8547 * <<
"----------------------------------------------------"
8554 * std::cerr << std::endl << std::endl
8555 * <<
"----------------------------------------------------"
8557 * std::cerr <<
"Unknown exception!" << std::endl
8558 * <<
"Aborting!" << std::endl
8559 * <<
"----------------------------------------------------"
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
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
cell_iterator end() const
unsigned int n_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
virtual types::global_cell_index n_global_active_cells() const override
types::subdomain_id locally_owned_subdomain() const override
virtual void execute_coarsening_and_refinement() override
virtual bool prepare_coarsening_and_refinement() override
virtual void clear() override
__global__ void set(Number *val, const Number s, const size_type N)
#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.
void reset_manifold(const types::manifold_id manifold_number)
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
@ 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.
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)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
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)
long double gamma(const unsigned int n)
int(& functions)(const void *v1, const void *v2)
const 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 > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
void swap(SmartPointer< T, P > &t1, SmartPointer< T, Q > &t2)
std::vector< unsigned int > vertices
types::boundary_id boundary_id
std::vector< std::vector< bool > > constant_modes
double aggregation_threshold
unsigned int smoother_sweeps
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
1.11.0