164 * <a name=
"Includefiles"></a>
165 * <h3>Include files</h3>
166 * The
set of include files is not much of a surprise any more at
this time:
169 * #include <deal.II/base/conditional_ostream.h>
170 * #include <deal.II/base/parameter_handler.h>
171 * #include <deal.II/base/utilities.h>
172 * #include <deal.II/base/index_set.h>
173 * #include <deal.II/base/quadrature_lib.h>
174 * #include <deal.II/base/function.h>
175 * #include <deal.II/base/logstream.h>
176 * #include <deal.II/base/timer.h>
177 * #include <deal.II/base/table_handler.h>
179 * #include <deal.II/lac/vector.h>
180 * #include <deal.II/lac/full_matrix.h>
181 * #include <deal.II/lac/sparsity_tools.h>
182 * #include <deal.II/lac/sparse_matrix.h>
183 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
184 * #include <deal.II/lac/block_sparsity_pattern.h>
185 * #include <deal.II/lac/solver_bicgstab.h>
186 * #include <deal.II/lac/precondition.h>
187 * #include <deal.II/lac/affine_constraints.h>
188 * #include <deal.II/lac/trilinos_sparse_matrix.h>
189 * #include <deal.II/lac/trilinos_block_sparse_matrix.h>
190 * #include <deal.II/lac/trilinos_vector.h>
191 * #include <deal.II/lac/trilinos_precondition.h>
192 * #include <deal.II/lac/trilinos_solver.h>
193 * #include <deal.II/lac/sparse_direct.h>
195 * #include <deal.II/grid/
tria.h>
196 * #include <deal.II/grid/grid_generator.h>
197 * #include <deal.II/grid/grid_refinement.h>
198 * #include <deal.II/grid/grid_tools.h>
199 * #include <deal.II/grid/tria_accessor.h>
200 * #include <deal.II/grid/tria_iterator.h>
201 * #include <deal.II/grid/grid_out.h>
202 * #include <deal.II/grid/manifold_lib.h>
204 * #include <deal.II/distributed/
tria.h>
205 * #include <deal.II/distributed/grid_refinement.h>
206 * #include <deal.II/distributed/solution_transfer.h>
208 * #include <deal.II/dofs/dof_handler.h>
209 * #include <deal.II/dofs/dof_accessor.h>
210 * #include <deal.II/dofs/dof_renumbering.h>
211 * #include <deal.II/dofs/dof_tools.h>
213 * #include <deal.II/fe/fe_q.h>
214 * #include <deal.II/fe/fe_system.h>
215 * #include <deal.II/fe/fe_values.h>
216 * #include <deal.II/fe/fe_dgq.h>
217 * #include <deal.II/fe/fe_tools.h>
219 * #include <deal.II/numerics/vector_tools.h>
220 * #include <deal.II/numerics/matrix_tools.h>
221 * #include <deal.II/numerics/data_out.h>
222 * #include <deal.II/numerics/error_estimator.h>
223 * #include <deal.II/numerics/fe_field_function.h>
224 * #include <deal.II/numerics/solution_transfer.h>
228 * And here the only two
new things among the header files: an include file in
229 * which
symmetric tensors of rank 2 and 4 are implemented, as introduced in
233 * #include <deal.II/base/symmetric_tensor.h>
237 * And a header that implements filters
for iterators looping over all
238 * cells. We will use
this when selecting only those cells
for output that are
239 * owned by the present process in a %
parallel program:
242 * #include <deal.II/grid/filtered_iterator.h>
245 * #include <iostream>
249 * This
final include file provides the <code>mkdir</code> function
250 * that we will use to create a directory
for output files,
if necessary:
253 * #include <sys/stat.h>
255 *
namespace ElastoPlastic
261 *
const unsigned int n_slices,
262 *
const double height,
267 *
Assert (input.n_levels() == 1,
268 * ExcMessage (
"The input triangulations must be coarse meshes."));
271 *
Assert(result.n_cells()==0, ExcMessage(
"resultin Triangulation need to be empty upon calling extrude_triangulation."));
272 *
Assert(height>0, ExcMessage(
"The height in extrude_triangulation needs to be positive."));
273 *
Assert(n_slices>=2, ExcMessage(
"The number of slices in extrude_triangulation needs to be at least 2."));
275 * std::vector<Point<3> > points(n_slices*input.n_used_vertices());
276 * std::vector<CellData<3> > cells;
277 * cells.reserve((n_slices-1)*input.n_active_cells());
279 *
for (
unsigned int slice=0; slice<n_slices; ++slice)
281 *
for (
unsigned int i=0; i<input.n_vertices(); ++i)
284 *
if (input.get_used_vertices()[i])
286 *
const Point<2> &v = input.get_vertices()[i];
287 * points[i+slice*input.n_vertices()](0) = v(0);
288 * points[i+slice*input.n_vertices()](1) = v(1);
289 * points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1);
295 * cell = input.begin_active(); cell != input.end(); ++cell)
297 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
300 *
for (
unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
303 * = cell->vertex_index(v)+slice*input.n_used_vertices();
305 * = cell->vertex_index(v)+(slice+1)*input.n_used_vertices();
308 * this_cell.material_id = cell->material_id();
309 * cells.push_back(this_cell);
315 * s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
317 * cell = input.begin_active(); cell != input.end(); ++cell)
320 *
for (
unsigned int f=0; f<4; ++f)
321 *
if (cell->at_boundary(f))
324 * bid =
std::max(bid, quad.boundary_id);
325 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
327 * quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_used_vertices();
328 * quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_used_vertices();
329 * quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_used_vertices();
330 * quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_used_vertices();
331 * s.boundary_quads.push_back(quad);
337 * cell = input.begin_active(); cell != input.end(); ++cell)
341 * quad.vertices[0] = cell->vertex_index(0);
342 * quad.vertices[1] = cell->vertex_index(1);
343 * quad.vertices[2] = cell->vertex_index(2);
344 * quad.vertices[3] = cell->vertex_index(3);
345 * s.boundary_quads.push_back(quad);
347 * quad.boundary_id = bid + 2;
348 *
for (
int i=0; i<4; ++i)
349 * quad.vertices[i] += (n_slices-1)*input.n_used_vertices();
350 * s.boundary_quads.push_back(quad);
353 * result.create_triangulation (points,
358 *
namespace Evaluation
370 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
371 * + stress[1][1]*stress[1][1]
372 * - stress[0][0]*stress[1][1]
373 * + 3*stress[0][1]*stress[0][1]);
374 * }
else if (dim == 3)
376 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
377 * + stress[1][1]*stress[1][1]
378 * + stress[2][2]*stress[2][2]
379 * - stress[0][0]*stress[1][1]
380 * - stress[1][1]*stress[2][2]
381 * - stress[0][0]*stress[2][2]
382 * + 3*( stress[0][1]*stress[0][1]
383 * +stress[1][2]*stress[1][2]
384 * +stress[0][2]*stress[0][2]) );
389 * -----------------------------------------------
390 *
"Perforated_strip_tension"
392 *
const double von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
393 * + stress[1][1]*stress[1][1]
394 * - stress[0][0]*stress[1][1]
395 * + 3*stress[0][1]*stress[0][1]);
396 * -----------------------------------------------
398 * plane strain / 3
d case
404 * -----------------------------------------------
416 *
return von_Mises_stress;
421 *
class PointValuesEvaluation
424 * PointValuesEvaluation (
const Point<dim> &evaluation_point);
432 * <<
"The evaluation point " << arg1
433 * <<
" was not found among the vertices of the present grid.");
440 * PointValuesEvaluation<dim>::
441 * PointValuesEvaluation (
const Point<dim> &evaluation_point)
443 * evaluation_point (evaluation_point)
450 * PointValuesEvaluation<dim>::
455 *
const unsigned int dofs_per_vertex = dof_handler.get_fe().dofs_per_vertex;
457 * ExcDimensionMismatch (
point_values.size(), dofs_per_vertex));
461 * cell = dof_handler.begin_active(),
462 * endc = dof_handler.end();
463 *
bool evaluation_point_found =
false;
464 *
for (; (cell!=endc) && !evaluation_point_found; ++cell)
466 *
if (cell->is_locally_owned() && !evaluation_point_found)
467 *
for (
unsigned int vertex=0;
468 * vertex<GeometryInfo<dim>::vertices_per_cell;
471 *
if (cell->vertex(vertex).distance (evaluation_point)
473 * cell->diameter() * 1
e-8)
475 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
477 *
point_values[id] = solution(cell->vertex_dof_index(vertex,
id));
480 * evaluation_point_found =
true;
487 * ExcEvaluationPointNotFound(evaluation_point));
496 * <a name=
"ThecodePointHistorycodeclass"></a>
497 * <h3>The <code>PointHistory</code>
class</h3>
501 * As was mentioned in the introduction, we have to store the old stress in
502 * quadrature
point so that we can compute the residual forces at
this point
503 * during the next time step. This alone would not warrant a structure with
504 * only one member, but in more complicated applications, we would have to
505 * store more information in quadrature points as well, such as the history
506 * variables of plasticity, etc. In essence, we have to store everything
507 * that affects the present state of the material here, which in plasticity
508 * is determined by the deformation history variables.
512 * We will not give
this class any meaningful functionality beyond being
513 * able to store data, i.e. there are no constructors, destructors, or other
514 * member
functions. In such cases of `dumb
' classes, we usually opt to
515 * declare them as <code>struct</code> rather than <code>class</code>, to
516 * indicate that they are closer to C-style structures than C++-style
521 * struct PointHistory
523 * SymmetricTensor<2,dim> old_stress;
524 * SymmetricTensor<2,dim> old_strain;
532 * <a name="ThecodeConstitutiveLawcodeclasstemplate"></a>
533 * <h3>The <code>ConstitutiveLaw</code> class template</h3>
537 * This class provides an interface for a constitutive law, i.e., for the
538 * relationship between strain @f$\varepsilon(\mathbf u)@f$ and stress
539 * @f$\sigma@f$. In this example we are using an elastoplastic material behavior
540 * with linear, isotropic hardening. Such materials are characterized by
541 * Young's modulus @f$E@f$, Poisson
's ratio @f$\nu@f$, the initial yield stress
542 * @f$\sigma_0@f$ and the isotropic hardening parameter @f$\gamma@f$. For @f$\gamma =
543 * 0@f$ we obtain perfect elastoplastic behavior.
547 * As explained in the paper that describes this program, the first Newton
548 * steps are solved with a completely elastic material model to avoid having
549 * to deal with both nonlinearities (plasticity and contact) at once. To this
550 * end, this class has a function <code>set_sigma_0()</code> that we use later
551 * on to simply set @f$\sigma_0@f$ to a very large value -- essentially
552 * guaranteeing that the actual stress will not exceed it, and thereby
553 * producing an elastic material. When we are ready to use a plastic model, we
554 * set @f$\sigma_0@f$ back to its proper value, using the same function. As a
555 * result of this approach, we need to leave <code>sigma_0</code> as the only
556 * non-const member variable of this class.
560 * class ConstitutiveLaw
563 * ConstitutiveLaw (const double E,
565 * const double sigma_0,
566 * const double gamma);
569 * set_sigma_0 (double sigma_zero);
572 * get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
573 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
576 * get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
577 * const std::vector<Tensor<2, dim> > &point_hessian,
578 * Tensor<5, dim> &stress_strain_tensor_grad) const;
581 * get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor,
582 * SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
583 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
586 * const double kappa;
589 * const double gamma;
591 * const SymmetricTensor<4, dim> stress_strain_tensor_kappa;
592 * const SymmetricTensor<4, dim> stress_strain_tensor_mu;
597 * The constructor of the ConstitutiveLaw class sets the required material
598 * parameter for our deformable body. Material parameters for elastic
599 * isotropic media can be defined in a variety of ways, such as the pair @f$E,
600 * \nu@f$ (elastic modulus and Poisson's number),
using the Lame parameters
601 * @f$\lambda,mu@f$ or several other commonly used conventions. Here, the
602 * constructor takes a description of material parameters in the form of
603 * @f$E,\nu@f$, but since
this turns out to these are not the coefficients that
604 * appear in the equations of the plastic projector, we immediately convert
605 * them into the more suitable set @f$\kappa,\mu@f$ of bulk and shear moduli. In
606 * addition, the constructor takes @f$\sigma_0@f$ (the yield stress absent any
607 * plastic strain) and @f$\gamma@f$ (the hardening parameter) as arguments. In
608 *
this constructor, we also compute the two principal components of the
609 * stress-strain relation and its linearization.
613 * ConstitutiveLaw<dim>::ConstitutiveLaw (
double E,
620 * --------------------
622 * kappa (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (3 * (1 - 2 * (nu / (1+nu))))),
623 * mu (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))),
624 * --------------------
625 * 3
d and plane strain
628 * kappa (E / (3 * (1 - 2 * nu))),
629 * mu (E / (2 * (1 + nu))),
632 * --------------------
637 * stress_strain_tensor_kappa (kappa
639 * unit_symmetric_tensor<dim>())),
640 * stress_strain_tensor_mu (2 * mu
641 * * (identity_tensor<dim>()
643 * unit_symmetric_tensor<dim>()) / 3.0))
649 * ConstitutiveLaw<dim>::set_sigma_0 (
double sigma_zero)
651 * sigma_0 = sigma_zero;
658 * <a name=
"ConstitutiveLawget_stress_strain_tensor"></a>
659 * <h4>ConstitutiveLaw::get_stress_strain_tensor</h4>
663 * This is the principal component of the constitutive law. It projects the
664 * deviatoric part of the stresses in a quadrature
point back to the yield
665 * stress (i.e., the original yield stress @f$\sigma_0@f$ plus the term that
666 * describes linear isotropic hardening). We need
this function to calculate
667 * the nonlinear residual in PlasticityContactProblem::residual_nl_system. The
668 * computations follow the formulas laid out in the introduction.
672 * The function returns whether the quadrature
point is plastic to allow
for
673 * some statistics
downstream on how many of the quadrature points are
674 * plastic and how many are elastic.
679 * ConstitutiveLaw<dim>::
684 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
690 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
693 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
695 * stress_strain_tensor = stress_strain_tensor_mu;
696 *
if (von_Mises_stress > sigma_0)
698 *
const double beta = sigma_0 / von_Mises_stress;
699 * stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
702 * stress_strain_tensor += stress_strain_tensor_kappa;
704 *
return (von_Mises_stress > sigma_0);
710 * ConstitutiveLaw<dim>::
716 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
719 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
721 *
if (von_Mises_stress > sigma_0)
724 *
const double deviator_strain_tensor_norm = deviator_strain_tensor.norm();
725 *
const double multiplier = -(1-
gamma)*sigma_0/(2*mu*
std::pow(deviator_strain_tensor_norm,3));
728 * multiplier_vector = 0;
730 *
for (
unsigned int i=0; i!=dim; ++i)
731 *
for (
unsigned int m=0; m!=dim; ++m)
732 *
for (
unsigned int n=0; n!=dim; ++n)
734 * multiplier_vector(i) += deviator_strain_tensor[m][n] *
735 * ( 0.5*( point_hessian[m][n][i] + point_hessian[n][m][i] )
736 * + ( m==n && dim==2 ? -1/dim*(point_hessian[0][0][i]
737 * + point_hessian[1][1][i]) : 0 )
738 * + ( m==n && dim==3 ? -1/dim*(point_hessian[0][0][i]
739 * + point_hessian[1][1][i]
740 * + point_hessian[2][2][i]) : 0 ) );
745 * -----------------------------------------------
746 *
"Perforated_strip_tension"
749 * -----------------------------------------------
751 * plane strain / 3
d case
754 *
const double VM_factor =
std::sqrt(1.5);
757 * -----------------------------------------------
763 *
for (
unsigned int i=0; i!=dim; ++i)
764 *
for (
unsigned int j=0; j!=dim; ++j)
765 *
for (
unsigned int k=0; k!=dim; ++k)
766 *
for (
unsigned int l=0;
l!=dim; ++
l)
767 *
for (
unsigned int m=0; m!=dim; ++m)
769 * stress_strain_tensor_grad[i][j][k][
l][m] = 1/VM_factor
771 * * stress_strain_tensor_mu[i][j][k][
l]
772 * * multiplier_vector(m);
778 * stress_strain_tensor_grad = 0;
781 *
return (von_Mises_stress > sigma_0);
788 * <a name=
"ConstitutiveLawget_linearized_stress_strain_tensors"></a>
789 * <h4>ConstitutiveLaw::get_linearized_stress_strain_tensors</h4>
793 * This function returns the linearized stress strain tensor, linearized
794 * around the solution @f$u^{i-1}@f$ of the previous Newton step @f$i-1@f$. The
795 * parameter <code>strain_tensor</code> (commonly denoted
796 * @f$\varepsilon(u^{i-1})@f$) must be passed as an argument, and serves as the
797 * linearization
point. The function returns the derivative of the nonlinear
798 * constitutive law in the variable stress_strain_tensor, as well as the
799 * stress-strain tensor of the linearized problem in
800 * stress_strain_tensor_linearized. See
801 * PlasticityContactProblem::assemble_nl_system where
this function is used.
806 * ConstitutiveLaw<dim>::
812 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
815 * stress_strain_tensor = stress_strain_tensor_mu;
816 * stress_strain_tensor_linearized = stress_strain_tensor_mu;
819 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
820 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
822 *
if (von_Mises_stress > sigma_0)
824 *
const double beta = sigma_0 / von_Mises_stress;
825 * stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
826 * stress_strain_tensor_linearized *= (
gamma + (1 -
gamma) * beta);
827 * deviator_stress_tensor /= deviator_stress_tensor_norm;
828 * stress_strain_tensor_linearized -= (1 -
gamma) * beta * 2 * mu
830 * deviator_stress_tensor);
833 * stress_strain_tensor += stress_strain_tensor_kappa;
834 * stress_strain_tensor_linearized += stress_strain_tensor_kappa;
839 * Finally, below we will need a function that computes the rotation
matrix
840 * induced by a displacement at a given
point. In fact, of course, the
841 * displacement at a single
point only has a direction and a magnitude, it
842 * is the change in direction and magnitude that induces rotations. In
844 * displacement, or, more specifically, from the curl.
848 * The formulas by which the rotation matrices are determined are a little
849 * awkward, especially in 3
d. For 2
d, there is a simpler way, so we
850 * implement
this function twice, once
for 2
d and once
for 3
d, so that we
851 * can compile and use the program in both space dimensions
if so desired --
852 * after all, deal.II is all about dimension independent programming and
853 * reuse of algorithm thoroughly tested with cheap computations in 2
d,
for
854 * the more expensive computations in 3
d. Here is one
case, where we have to
855 * implement different algorithms
for 2
d and 3
d, but then can write the rest
856 * of the program in a way that is independent of the space dimension.
860 * So, without further ado to the 2
d implementation:
864 * get_rotation_matrix (
const std::vector<
Tensor<1,2> > &grad_u)
868 * First, compute the curl of the velocity field from the
gradients. Note
869 * that we are in 2
d, so the rotation is a
scalar:
872 *
const double curl = (grad_u[1][0] - grad_u[0][1]);
876 * From
this, compute the
angle of rotation:
879 *
const double angle = std::atan (curl);
883 * And from
this, build the antisymmetric rotation
matrix:
886 *
const double t[2][2] = {{
cos(angle),
sin(angle) },
887 * {-
sin(angle),
cos(angle) }
895 * The 3
d case is a little more contrived:
899 * get_rotation_matrix (
const std::vector<
Tensor<1,3> > &grad_u)
903 * Again
first compute the curl of the velocity field. This time, it is a
907 *
const Point<3> curl (grad_u[2][1] - grad_u[1][2],
908 * grad_u[0][2] - grad_u[2][0],
909 * grad_u[1][0] - grad_u[0][1]);
913 * From
this vector,
using its magnitude, compute the tangent of the
angle
914 * of rotation, and from it the actual
angle:
917 *
const double tan_angle =
std::sqrt(curl*curl);
918 *
const double angle = std::atan (tan_angle);
922 * Now, here
's one problem: if the angle of rotation is too small, that
923 * means that there is no rotation going on (for example a translational
924 * motion). In that case, the rotation matrix is the identity matrix.
928 * The reason why we stress that is that in this case we have that
929 * <code>tan_angle==0</code>. Further down, we need to divide by that
930 * number in the computation of the axis of rotation, and we would get
931 * into trouble when dividing doing so. Therefore, let's shortcut
this and
938 *
static const double rotation[3][3]
939 * = {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } };
946 * Otherwise compute the real rotation
matrix. The algorithm
for this is
947 * not exactly obvious, but can be found in a number of books,
948 * particularly on computer games where rotation is a very frequent
949 * operation. Online, you can find a description at
951 * the signs as here) at
957 *
const double t = 1-c;
959 *
const Point<3> axis = curl/tan_angle;
960 *
const double rotation[3][3]
962 * t *axis[0] *axis[0]+c,
963 * t *axis[0] *axis[1]+s *axis[2],
964 * t *axis[0] *axis[2]-s *axis[1]
967 * t *axis[0] *axis[1]-s *axis[2],
968 * t *axis[1] *axis[1]+c,
969 * t *axis[1] *axis[2]+s *axis[0]
972 * t *axis[0] *axis[2]+s *axis[1],
973 * t *axis[1] *axis[1]-s *axis[0],
974 * t *axis[2] *axis[2]+c
983 * <h3>Equation data: Body forces, boundary forces,
984 * incremental boundary
values</h3>
988 * The following should be relatively standard. We need classes
for
989 * the boundary forcing term (which we here choose to be zero)
990 * and incremental boundary
values.
993 *
namespace EquationData
1239 * ----------------------------- TimoshenkoBeam ---------------------------------------
1442 * ------------------------- Thick_tube_internal_pressure ----------------------------------
1611 * ------------------------- Perforated_strip_tension ----------------------------------
1781 * ------------------------- Cantiliver_beam_3d ----------------------------------
1784 *
template <
int dim>
1785 *
class IncrementalBoundaryForce :
public Function<dim>
1788 * IncrementalBoundaryForce (
const double present_time,
1789 *
const double end_time);
1801 *
const double present_time,
1807 *
template <
int dim>
1808 * IncrementalBoundaryForce<dim>::
1809 * IncrementalBoundaryForce (
const double present_time,
1810 *
const double end_time)
1813 * present_time (present_time),
1814 * end_time (end_time),
1819 *
template <
int dim>
1821 * IncrementalBoundaryForce<dim>::vector_value (
const Point<dim> &p,
1826 * ExcDimensionMismatch (
values.size(), dim));
1828 *
const double eps = 1.e-7 * height;
1832 * pressure should be imposed on the top surface, y = height
1843 * compute the fraction of imposed force
1846 *
const double frac = present_time/end_time;
1851 *
template <
int dim>
1853 * IncrementalBoundaryForce<dim>::
1854 * vector_value_list (
const std::vector<
Point<dim> > &points,
1857 *
const unsigned int n_points = points.size();
1859 *
Assert (value_list.size() == n_points,
1860 * ExcDimensionMismatch (value_list.size(), n_points));
1862 *
for (
unsigned int p=0; p<n_points; ++p)
1863 * IncrementalBoundaryForce<dim>::vector_value (points[p], value_list[p]);
1867 *
template <
int dim>
1871 * BodyForce () :
Functions::ZeroFunction<dim> (dim) {}
1875 *
template <
int dim>
1876 *
class IncrementalBoundaryValues :
public Function<dim>
1879 * IncrementalBoundaryValues (
const double present_time,
1880 *
const double end_time);
1893 *
const double present_time,
1898 *
template <
int dim>
1899 * IncrementalBoundaryValues<dim>::
1900 * IncrementalBoundaryValues (
const double present_time,
1901 *
const double end_time)
1904 * present_time (present_time),
1905 * end_time (end_time)
1909 *
template <
int dim>
1911 * IncrementalBoundaryValues<dim>::
1916 * ExcDimensionMismatch (
values.size(), dim));
1923 *
template <
int dim>
1925 * IncrementalBoundaryValues<dim>::
1926 * vector_value_list (
const std::vector<
Point<dim> > &points,
1929 *
const unsigned int n_points = points.size();
1931 *
Assert (value_list.size() == n_points,
1932 * ExcDimensionMismatch (value_list.size(), n_points));
1934 *
for (
unsigned int p=0; p<n_points; ++p)
1935 * IncrementalBoundaryValues<dim>::vector_value (points[p], value_list[p]);
1940 * -------------------------------------------------------------------------------
1946 *
namespace DualFunctional
1949 *
template <
int dim>
1957 *
const ConstitutiveLaw<dim> &constitutive_law,
1963 *
template <
int dim>
1964 *
class PointValuesEvaluation :
public DualFunctionalBase<dim>
1967 * PointValuesEvaluation (
const Point<dim> &evaluation_point);
1973 *
const ConstitutiveLaw<dim> &constitutive_law,
1979 * <<
"The evaluation point " << arg1
1980 * <<
" was not found among the vertices of the present grid.");
1987 *
template <
int dim>
1988 * PointValuesEvaluation<dim>::
1989 * PointValuesEvaluation (
const Point<dim> &evaluation_point)
1991 * evaluation_point (evaluation_point)
1995 *
template <
int dim>
1997 * PointValuesEvaluation<dim>::
2000 *
const ConstitutiveLaw<dim> &,
2004 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2005 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2008 * cell_dual = dof_handler_dual.begin_active(),
2009 * endc_dual = dof_handler_dual.end();
2010 *
for (; cell_dual!=endc_dual; ++cell_dual)
2011 *
for (
unsigned int vertex=0;
2012 * vertex<GeometryInfo<dim>::vertices_per_cell;
2014 *
if (cell_dual->vertex(vertex).distance(evaluation_point)
2017 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
2019 * rhs_dual(cell_dual->vertex_dof_index(vertex,
id)) = 1;
2024 *
AssertThrow (
false, ExcEvaluationPointNotFound(evaluation_point));
2028 *
template <
int dim>
2029 *
class PointXDerivativesEvaluation :
public DualFunctionalBase<dim>
2032 * PointXDerivativesEvaluation (
const Point<dim> &evaluation_point);
2038 *
const ConstitutiveLaw<dim> &constitutive_law,
2044 * <<
"The evaluation point " << arg1
2045 * <<
" was not found among the vertices of the present grid.");
2052 *
template <
int dim>
2053 * PointXDerivativesEvaluation<dim>::
2054 * PointXDerivativesEvaluation (
const Point<dim> &evaluation_point)
2056 * evaluation_point (evaluation_point)
2060 *
template <
int dim>
2062 * PointXDerivativesEvaluation<dim>::
2065 *
const ConstitutiveLaw<dim> &,
2069 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2072 *
FEValues<dim> fe_values (dof_handler_dual.get_fe(), quadrature,
2076 *
const unsigned int n_q_points = fe_values.n_quadrature_points;
2077 *
Assert ( n_q_points==quadrature.size() , ExcInternalError() );
2078 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2081 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2083 *
double total_volume = 0;
2086 * cell = dof_handler_dual.begin_active(),
2087 * endc = dof_handler_dual.end();
2088 *
for (; cell!=endc; ++cell)
2089 *
if (cell->center().distance(evaluation_point) <=
2092 * fe_values.reinit (cell);
2095 *
for (
unsigned int q=0; q<n_q_points; ++q)
2097 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2099 * cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
2100 * fe_values.JxW (q);
2103 * total_volume += fe_values.JxW (q);
2106 * cell->get_dof_indices (local_dof_indices);
2107 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2109 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2114 * ExcEvaluationPointNotFound(evaluation_point));
2116 * rhs_dual *= 1./total_volume;
2121 *
template <
int dim>
2122 *
class MeanDisplacementFace :
public DualFunctionalBase<dim>
2125 * MeanDisplacementFace (
const unsigned int face_id,
2126 *
const std::vector<bool> comp_mask);
2132 *
const ConstitutiveLaw<dim> &constitutive_law,
2137 *
const unsigned int face_id;
2138 *
const std::vector<bool> comp_mask;
2142 *
template <
int dim>
2143 * MeanDisplacementFace<dim>::
2144 * MeanDisplacementFace (
const unsigned int face_id,
2145 *
const std::vector<bool> comp_mask )
2147 * face_id (face_id),
2148 * comp_mask (comp_mask)
2151 * ExcDimensionMismatch (comp_mask.size(), dim) );
2155 *
template <
int dim>
2157 * MeanDisplacementFace<dim>::
2160 *
const ConstitutiveLaw<dim> &,
2166 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2168 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2172 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2173 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2174 *
const unsigned int n_face_q_points = face_quadrature.size();
2177 * ExcDimensionMismatch (dofs_per_vertex, dim) );
2179 * std::vector<unsigned int> comp_vector(dofs_per_vertex);
2180 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2184 * comp_vector[i] = 1;
2190 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2194 * bound_size : size of the boundary, in 2
d is the length
2195 * and in the 3
d case, area
2198 *
double bound_size = 0.;
2201 * cell = dof_handler_dual.begin_active(),
2202 * endc = dof_handler_dual.end();
2203 *
bool evaluation_face_found =
false;
2204 *
for (; cell!=endc; ++cell)
2207 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2209 *
if (cell->face(face)->at_boundary()
2211 * cell->face(face)->boundary_id() == face_id)
2213 *
if (!evaluation_face_found)
2215 * evaluation_face_found =
true;
2217 * fe_face_values.reinit (cell, face);
2219 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2221 * bound_size += fe_face_values.JxW(q_point);
2223 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2225 *
const unsigned int
2226 * component_i = dof_handler_dual.get_fe().system_to_component_index(i).first;
2228 * cell_rhs(i) += (fe_face_values.shape_value(i,q_point) *
2229 * comp_vector[component_i] *
2230 * fe_face_values.JxW(q_point));
2238 * cell->get_dof_indices (local_dof_indices);
2239 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2241 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2246 *
AssertThrow(evaluation_face_found, ExcInternalError());
2248 * rhs_dual /= bound_size;
2253 *
template <
int dim>
2254 *
class MeanStressFace :
public DualFunctionalBase<dim>
2257 * MeanStressFace (
const unsigned int face_id,
2258 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2264 *
const ConstitutiveLaw<dim> &constitutive_law,
2269 *
const unsigned int face_id;
2270 *
const std::vector<std::vector<unsigned int> > comp_stress;
2274 *
template <
int dim>
2275 * MeanStressFace<dim>::
2276 * MeanStressFace (
const unsigned int face_id,
2277 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2279 * face_id (face_id),
2280 * comp_stress (comp_stress)
2283 * ExcDimensionMismatch (comp_stress.size(), dim) );
2287 *
template <
int dim>
2289 * MeanStressFace<dim>::
2292 *
const ConstitutiveLaw<dim> &constitutive_law,
2298 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2300 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2304 *
FEFaceValues<dim> fe_face_values_dual (dof_handler_dual.get_fe(), face_quadrature,
2307 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2308 *
const unsigned int n_face_q_points = face_quadrature.size();
2310 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_face_q_points);
2315 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2319 * bound_size : size of the boundary, in 2
d is the length
2320 * and in the 3
d case, area
2323 *
double bound_size = 0.;
2325 *
bool evaluation_face_found =
false;
2328 * cell_dual = dof_handler_dual.begin_active(),
2329 * endc_dual = dof_handler_dual.end(),
2330 * cell = dof_handler.begin_active();
2334 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2337 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2339 *
if (cell_dual->face(face)->at_boundary()
2341 * cell_dual->face(face)->boundary_id() == face_id)
2343 *
if (!evaluation_face_found)
2345 * evaluation_face_found =
true;
2348 * fe_face_values.reinit (cell, face);
2349 * fe_face_values_dual.reinit (cell_dual, face);
2351 * fe_face_values[displacement].get_function_symmetric_gradients(solution,
2354 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2356 * bound_size += fe_face_values_dual.JxW(q_point);
2358 * constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2359 * stress_strain_tensor);
2361 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2364 * stress_phi_i = stress_strain_tensor
2365 * * fe_face_values_dual[displacement].symmetric_gradient(i, q_point);
2367 *
for (
unsigned int k=0; k!=dim; ++k)
2369 *
for (
unsigned int l=0;
l!=dim; ++
l)
2371 *
if ( comp_stress[k][l] == 1 )
2373 * cell_rhs(i) += stress_phi_i[k][
l]
2375 * fe_face_values_dual.JxW(q_point);
2388 * cell_dual->get_dof_indices (local_dof_indices);
2389 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2391 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2396 *
AssertThrow(evaluation_face_found, ExcInternalError());
2398 * rhs_dual /= bound_size;
2403 *
template <
int dim>
2404 *
class MeanStressDomain :
public DualFunctionalBase<dim>
2407 * MeanStressDomain (
const std::string &base_mesh,
2408 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2414 *
const ConstitutiveLaw<dim> &constitutive_law,
2419 *
const std::string base_mesh;
2420 *
const std::vector<std::vector<unsigned int> > comp_stress;
2424 *
template <
int dim>
2425 * MeanStressDomain<dim>::
2426 * MeanStressDomain (
const std::string &base_mesh,
2427 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2429 * base_mesh (base_mesh),
2430 * comp_stress (comp_stress)
2433 * ExcDimensionMismatch (comp_stress.size(), dim) );
2437 *
template <
int dim>
2439 * MeanStressDomain<dim>::
2442 *
const ConstitutiveLaw<dim> &constitutive_law,
2446 *
AssertThrow (base_mesh ==
"Cantiliver_beam_3d", ExcNotImplemented());
2451 * Mean stress at the specified domain is of interest.
2452 * The interest domains are located on the bottom and top of the flanges
2453 * close to the clamped face, z = 0
2454 * top domain: height/2 - thickness_flange <= y <= height/2
2455 * 0 <= z <= 2 * thickness_flange
2456 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
2457 * 0 <= z <= 2 * thickness_flange
2463 *
const double height = 200
e-3,
2464 * thickness_flange = 10
e-3;
2466 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2468 *
const QGauss<dim> quadrature_formula(dof_handler_dual.get_fe().tensor_degree()+1);
2470 *
FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
2472 *
FEValues<dim> fe_values_dual (dof_handler_dual.get_fe(), quadrature_formula,
2475 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2476 *
const unsigned int n_q_points = quadrature_formula.size();
2478 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
2483 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2487 * domain_size : size of the interested domain, in 2
d is the area
2491 *
double domain_size = 0.;
2493 *
bool evaluation_domain_found =
false;
2496 * cell_dual = dof_handler_dual.begin_active(),
2497 * endc_dual = dof_handler_dual.end(),
2498 * cell = dof_handler.begin_active();
2502 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2504 *
const double y = cell->center()[1],
2505 * z = cell->center()[2];
2508 * top domain: height/2 - thickness_flange <= y <= height/2
2509 * 0 <= z <= 2 * thickness_flange
2510 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
2511 * 0 <= z <= 2 * thickness_flange
2514 *
if ( ((z > 0) && (z < 2*thickness_flange)) &&
2515 * ( ((y > height/2 - thickness_flange) && (y < height/2)) ||
2516 * ((y > -height/2) && (y < -height/2 + thickness_flange)) ) )
2520 *
if (!evaluation_domain_found)
2522 * evaluation_domain_found =
true;
2525 * fe_values.reinit(cell);
2526 * fe_values_dual.reinit(cell_dual);
2528 * fe_values[displacement].get_function_symmetric_gradients(solution,
2531 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
2533 * domain_size += fe_values_dual.JxW(q_point);
2535 * constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2536 * stress_strain_tensor);
2538 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2541 * stress_phi_i = stress_strain_tensor
2542 * * fe_values_dual[displacement].symmetric_gradient(i, q_point);
2544 *
for (
unsigned int k=0; k!=dim; ++k)
2546 *
for (
unsigned int l=0;
l!=dim; ++
l)
2548 *
if ( comp_stress[k][l] == 1 )
2550 * cell_rhs(i) += stress_phi_i[k][
l]
2552 * fe_values_dual.JxW(q_point);
2564 * cell_dual->get_dof_indices (local_dof_indices);
2565 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2567 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2572 *
AssertThrow(evaluation_domain_found, ExcInternalError());
2574 * rhs_dual /= domain_size;
2579 *
template <
int dim>
2580 *
class MeanStrainEnergyFace :
public DualFunctionalBase<dim>
2583 * MeanStrainEnergyFace (
const unsigned int face_id,
2587 *
void assemble_rhs_nonlinear (
const DoFHandler<dim> &primal_dof_handler,
2593 *
const unsigned int face_id;
2599 *
template <
int dim>
2600 * MeanStrainEnergyFace<dim>::
2601 * MeanStrainEnergyFace (
const unsigned int face_id,
2605 * face_id (face_id),
2606 * lambda_function (&lambda_function),
2607 * mu_function (&mu_function)
2611 *
template <
int dim>
2613 * MeanStrainEnergyFace<dim>::
2621 * Assemble right hand side of the dual problem when the quantity of interest is
2622 * a nonlinear functinoal. In
this case, the QoI should be linearized which depends
2623 * on the solution of the primal problem.
2624 * The extracter of the linearized QoI functional is the
gradient of the the original
2625 * QoI functional with the primal solution
values.
2633 * rhs.reinit (dof_handler.n_dofs());
2635 *
const QGauss<dim-1> face_quadrature(dof_handler.get_fe().tensor_degree()+1);
2636 *
FEFaceValues<dim> primal_fe_face_values (primal_dof_handler.get_fe(), face_quadrature,
2644 *
const unsigned int dofs_per_vertex = primal_dof_handler.get_fe().dofs_per_vertex;
2645 *
const unsigned int n_face_q_points = face_quadrature.size();
2646 *
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
2649 * ExcDimensionMismatch (dofs_per_vertex, dim) );
2651 * std::vector< std::vector< Tensor<1,dim> > > primal_solution_gradients;
2652 * primal_solution_gradients.resize(n_face_q_points);
2654 * std::vector<std::vector<Tensor<2,dim> > > primal_solution_hessians;
2655 * primal_solution_hessians.resize (n_face_q_points);
2657 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
2659 * primal_solution_gradients[i].resize (dofs_per_vertex);
2660 * primal_solution_hessians[i].resize (dofs_per_vertex);
2663 * std::vector<double> lambda_values (n_face_q_points);
2664 * std::vector<double> mu_values (n_face_q_points);
2668 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2672 * bound_size : size of the boundary, in 2
d is the length
2673 * and in the 3
d case, area
2676 *
double bound_size = 0.;
2678 *
bool evaluation_face_found =
false;
2681 * primal_cell = primal_dof_handler.begin_active(),
2682 * primal_endc = primal_dof_handler.end();
2685 * cell = dof_handler.begin_active(),
2686 * endc = dof_handler.end();
2688 *
for (; cell!=endc; ++cell, ++primal_cell)
2691 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2693 *
if (cell->face(face)->at_boundary()
2695 * cell->face(face)->boundary_id() == face_id)
2697 *
if (!evaluation_face_found)
2699 * evaluation_face_found =
true;
2701 * primal_fe_face_values.reinit (primal_cell, face);
2703 * primal_fe_face_values.get_function_gradients (primal_solution,
2704 * primal_solution_gradients);
2706 * primal_fe_face_values.get_function_hessians (primal_solution,
2707 * primal_solution_hessians);
2709 * lambda_function->value_list (primal_fe_face_values.get_quadrature_points(), lambda_values);
2710 * mu_function->value_list (primal_fe_face_values.get_quadrature_points(), mu_values);
2712 * fe_face_values.reinit (cell, face);
2714 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2716 * bound_size += primal_fe_face_values.JxW(q_point);
2718 *
for (
unsigned int m=0; m<dofs_per_cell; ++m)
2720 *
const unsigned int
2721 * component_m = dof_handler.get_fe().system_to_component_index(m).first;
2723 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2725 *
for (
unsigned int j=0; j!=dofs_per_vertex; ++j)
2727 * cell_rhs(m) += fe_face_values.shape_value(m,q_point) *
2729 * lambda_values[q_point] *
2731 * primal_solution_hessians[q_point][i][i][component_m] * primal_solution_gradients[q_point][j][j]
2733 * primal_solution_gradients[q_point][i][i] * primal_solution_hessians[q_point][j][j][component_m]
2736 * mu_values[q_point] *
2738 * 2*primal_solution_hessians[q_point][j][i][component_m] * primal_solution_gradients[q_point][j][i]
2740 * primal_solution_hessians[q_point][i][j][component_m] * primal_solution_gradients[q_point][j][i]
2742 * primal_solution_gradients[q_point][i][j] * primal_solution_hessians[q_point][j][i][component_m]
2745 * primal_fe_face_values.JxW(q_point);
2758 * cell->get_dof_indices (local_dof_indices);
2759 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2761 * rhs(local_dof_indices[i]) += cell_rhs(i);
2766 *
AssertThrow(evaluation_face_found, ExcInternalError());
2768 * rhs *= 1./(2*bound_size);
2781 *
template <
int dim>
2788 *
const ConstitutiveLaw<dim> &constitutive_law,
2789 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2790 *
const unsigned int ×tep_no,
2791 *
const std::string &output_dir,
2792 *
const std::string &base_mesh,
2793 *
const double &present_time,
2794 *
const double &end_time);
2796 *
void compute_error_DWR (
Vector<float> &estimated_error_per_cell);
2801 *
void setup_system ();
2802 *
void compute_dirichlet_constraints ();
2803 *
void assemble_matrix ();
2804 *
void assemble_rhs ();
2806 *
void output_results ();
2812 *
const unsigned int fe_degree;
2815 *
const unsigned int fe_degree_dual;
2820 *
const QGauss<dim - 1> face_quadrature_formula;
2830 *
const ConstitutiveLaw<dim> constitutive_law;
2835 *
unsigned int timestep_no;
2836 * std::string output_dir;
2837 *
const std::string base_mesh;
2838 *
double present_time;
2848 *
const ConstitutiveLaw<dim> &constitutive_law,
2849 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2850 *
const unsigned int ×tep_no,
2851 *
const std::string &output_dir,
2852 *
const std::string &base_mesh,
2853 *
const double &present_time,
2854 *
const double &end_time)
2858 * solution(solution),
2859 * fe_degree(fe.tensor_degree()),
2860 * fe_degree_dual(fe_degree + 1),
2861 * fe_dual(
FE_Q<dim>(fe_degree_dual), dim),
2863 * quadrature_formula (fe_degree_dual + 1),
2864 * face_quadrature_formula (fe_degree_dual + 1),
2865 * constitutive_law (constitutive_law),
2867 * dual_functional (&dual_functional),
2868 * timestep_no (timestep_no),
2869 * output_dir (output_dir),
2870 * base_mesh (base_mesh),
2871 * present_time (present_time),
2872 * end_time (end_time)
2877 * DualSolver<dim>::~DualSolver()
2879 * dof_handler_dual.clear ();
2884 *
void DualSolver<dim>::setup_system()
2886 * dof_handler.distribute_dofs(fe);
2888 * dof_handler_dual.distribute_dofs (fe_dual);
2889 * std::cout <<
" Number of degrees of freedom in dual problem: "
2890 * << dof_handler_dual.n_dofs()
2893 * constraints_hanging_nodes_dual.clear ();
2895 * constraints_hanging_nodes_dual);
2896 * constraints_hanging_nodes_dual.close ();
2898 * compute_dirichlet_constraints();
2900 * sparsity_pattern_dual.reinit (dof_handler_dual.n_dofs(),
2901 * dof_handler_dual.n_dofs(),
2902 * dof_handler_dual.max_couplings_between_dofs());
2907 * constraints_hanging_nodes_dual.condense (sparsity_pattern_dual);
2910 * constraints_dirichlet_and_hanging_nodes_dual.condense (sparsity_pattern_dual);
2912 * sparsity_pattern_dual.compress();
2914 * system_matrix_dual.reinit (sparsity_pattern_dual);
2916 * solution_dual.reinit (dof_handler_dual.n_dofs());
2917 * system_rhs_dual.reinit (dof_handler_dual.n_dofs());
2922 *
void DualSolver<dim>::compute_dirichlet_constraints()
2924 * constraints_dirichlet_and_hanging_nodes_dual.clear ();
2925 * constraints_dirichlet_and_hanging_nodes_dual.merge(constraints_hanging_nodes_dual);
2927 * std::vector<bool> component_mask(dim);
2929 *
if (base_mesh ==
"Timoshenko beam")
2933 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2934 * constraints_dirichlet_and_hanging_nodes_dual,
2937 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
2941 * the boundary x = 0
2944 * component_mask[0] =
true;
2945 * component_mask[1] =
false;
2948 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2949 * constraints_dirichlet_and_hanging_nodes_dual,
2953 * the boundary y = 0
2956 * component_mask[0] =
false;
2957 * component_mask[1] =
true;
2960 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2961 * constraints_dirichlet_and_hanging_nodes_dual,
2964 *
else if (base_mesh ==
"Perforated_strip_tension")
2968 * the boundary x = 0
2971 * component_mask[0] =
true;
2972 * component_mask[1] =
false;
2973 * component_mask[2] =
false;
2976 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2977 * constraints_dirichlet_and_hanging_nodes_dual,
2981 * the boundary y = 0
2984 * component_mask[0] =
false;
2985 * component_mask[1] =
true;
2986 * component_mask[2] =
false;
2989 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2990 * constraints_dirichlet_and_hanging_nodes_dual,
2994 * the boundary y = imposed incremental displacement
2997 * component_mask[0] =
false;
2998 * component_mask[1] =
true;
2999 * component_mask[2] =
false;
3002 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
3003 * constraints_dirichlet_and_hanging_nodes_dual,
3006 *
else if (base_mesh ==
"Cantiliver_beam_3d")
3010 * the boundary x = y = z = 0
3013 * component_mask[0] =
true;
3014 * component_mask[1] =
true;
3015 * component_mask[2] =
true;
3018 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
3019 * constraints_dirichlet_and_hanging_nodes_dual,
3027 * constraints_dirichlet_and_hanging_nodes_dual.close();
3032 *
void DualSolver<dim>::assemble_matrix()
3039 *
const unsigned int dofs_per_cell_dual = fe_dual.dofs_per_cell;
3040 *
const unsigned int n_q_points = quadrature_formula.size();
3044 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell_dual);
3047 * cell_dual = dof_handler_dual.begin_active(),
3048 * endc_dual = dof_handler_dual.end(),
3049 * cell = dof_handler.begin_active();
3053 *
for (; cell_dual != endc_dual; ++cell_dual, ++cell)
3054 *
if (cell_dual->is_locally_owned())
3056 * fe_values.reinit(cell);
3058 * fe_values_dual.reinit(cell_dual);
3061 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3062 * fe_values[displacement].get_function_symmetric_gradients(solution,
3065 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3069 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3070 * stress_strain_tensor_linearized,
3071 * stress_strain_tensor);
3073 *
for (
unsigned int i = 0; i < dofs_per_cell_dual; ++i)
3076 * stress_phi_i = stress_strain_tensor_linearized
3077 * * fe_values_dual[displacement].symmetric_gradient(i, q_point);
3079 *
for (
unsigned int j = 0; j < dofs_per_cell_dual; ++j)
3081 * * fe_values_dual[displacement].symmetric_gradient(j, q_point)
3082 * * fe_values_dual.JxW(q_point));
3088 * cell_dual->get_dof_indices(local_dof_indices);
3089 * constraints_dirichlet_and_hanging_nodes_dual.distribute_local_to_global(cell_matrix,
3090 * local_dof_indices,
3091 * system_matrix_dual);
3099 *
void DualSolver<dim>::assemble_rhs()
3101 * dual_functional->assemble_rhs (dof_handler, solution, constitutive_law,
3102 * dof_handler_dual, system_rhs_dual);
3103 * constraints_dirichlet_and_hanging_nodes_dual.condense (system_rhs_dual);
3108 *
void DualSolver<dim>::solve()
3112 * +++ direct solver +++++++++
3120 * After the decomposition, we can use A_direct like a
matrix representing
3121 * the inverse of our system
matrix, so to compute the solution we just
3122 * have to multiply with the right hand side vector:
3125 * A_direct.vmult(solution_dual, system_rhs_dual);
3129 * ++++ iterative solver ++ CG ++++ doesn
't work
3130 * SolverControl solver_control (5000, 1e-12);
3131 * SolverCG<> cg (solver_control);
3135 * PreconditionSSOR<> preconditioner;
3136 * preconditioner.initialize(system_matrix_dual, 1.2);
3140 * cg.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3145 * ++++ iterative solver ++ BiCGStab ++++++ doesn't work
3152 * preconditioner.
initialize(system_matrix_dual, 1.0);
3156 * bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3161 * +++++++++++++++++++++++++++++++++++++++++++++++++
3167 * constraints_dirichlet_and_hanging_nodes_dual.distribute (solution_dual);
3171 *
void DualSolver<dim>::output_results()
3173 * std::string filename = (output_dir +
"dual-solution-" +
3175 * std::ofstream output (filename.c_str());
3178 * std::vector<std::string> solution_names;
3182 * solution_names.push_back (
"displacement");
3185 * solution_names.push_back (
"x_displacement");
3186 * solution_names.push_back (
"y_displacement");
3189 * solution_names.push_back (
"x_displacement");
3190 * solution_names.push_back (
"y_displacement");
3191 * solution_names.push_back (
"z_displacement");
3194 *
Assert (
false, ExcNotImplemented());
3196 * data_out.add_data_vector (solution_dual, solution_names);
3197 * data_out.build_patches ();
3198 * data_out.write_vtk (output);
3202 *
void DualSolver<dim>::compute_error_DWR (
Vector<float> &estimated_error_per_cell)
3205 * ExcDimensionMismatch (estimated_error_per_cell.size(),
triangulation->n_global_active_cells()));
3209 * solve the dual problem
3213 * assemble_matrix ();
3216 * output_results ();
3220 * compuate the dual weights
3227 * constraints_dirichlet_and_hanging_nodes_dual,
3232 * constraints_hanging_nodes);
3233 * constraints_hanging_nodes.close();
3236 * constraints_dirichlet_and_hanging_nodes_dual,
3239 * constraints_hanging_nodes,
3244 * estimate the error
3254 *
const unsigned int n_q_points = quadrature_formula.size();
3255 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3259 * std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
3260 *
for (
unsigned int i=0; i!=n_q_points; ++i)
3262 * cell_hessians[i].resize (dim);
3264 * std::vector<Vector<double> > dual_weights_cell_values (n_q_points,
Vector<double>(dim));
3266 *
const EquationData::BodyForce<dim> body_force;
3267 * std::vector<Vector<double> > body_force_values (n_q_points,
Vector<double>(dim));
3277 * fe_face_values_neighbor (fe_dual, face_quadrature_formula,
3285 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
3286 * std::vector<Vector<double> > jump_residual (n_face_q_points,
Vector<double>(dim));
3287 * std::vector<Vector<double> > dual_weights_face_values (n_face_q_points,
Vector<double>(dim));
3289 * std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
3290 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3292 * cell_grads[i].resize (dim);
3294 * std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
3295 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3297 * neighbor_grads[i].resize (dim);
3305 *
typename std::map<typename DoFHandler<dim>::face_iterator,
Vector<double> >
3308 * cell = dof_handler_dual.begin_active(),
3309 * endc = dof_handler_dual.
end();
3310 *
for (; cell!=endc; ++cell)
3311 *
if (cell->is_locally_owned())
3313 *
for (
unsigned int face_no=0;
3314 * face_no<GeometryInfo<dim>::faces_per_cell;
3317 * face_integrals[cell->face(face_no)].reinit (dim);
3318 * face_integrals[cell->face(face_no)] = -1e20;
3322 * std::vector<Vector<float> > error_indicators_vector;
3323 * error_indicators_vector.resize(
triangulation->n_active_cells(),
3328 * ----------------- estimate_some -------------------------
3331 * cell = dof_handler_dual.begin_active();
3332 *
unsigned int present_cell = 0;
3333 *
for (; cell!=endc; ++cell, ++present_cell)
3334 *
if (cell->is_locally_owned())
3338 * --------------- integrate_over_cell -------------------
3341 * fe_values.reinit(cell);
3342 * body_force.vector_value_list(fe_values.get_quadrature_points(),
3343 * body_force_values);
3344 * fe_values[displacement].get_function_symmetric_gradients(primal_solution,
3346 * fe_values.get_function_hessians(primal_solution, cell_hessians);
3348 * fe_values.get_function_values(dual_weights,
3349 * dual_weights_cell_values);
3351 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3353 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3354 * stress_strain_tensor_linearized,
3355 * stress_strain_tensor);
3356 * constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
3357 * cell_hessians[q_point],
3358 * stress_strain_tensor_grad);
3360 *
for (
unsigned int i=0; i!=dim; ++i)
3362 * error_indicators_vector[present_cell](i) +=
3363 * body_force_values[q_point](i)*
3364 * dual_weights_cell_values[q_point](i)*
3365 * fe_values.JxW(q_point);
3366 *
for (
unsigned int j=0; j!=dim; ++j)
3368 *
for (
unsigned int k=0; k!=dim; ++k)
3370 *
for (
unsigned int l=0;
l!=dim; ++
l)
3372 * error_indicators_vector[present_cell](i) +=
3373 * ( stress_strain_tensor[i][j][k][l]*
3374 * 0.5*(cell_hessians[q_point][k][l][j]
3376 * cell_hessians[q_point][l][k][j])
3377 * + stress_strain_tensor_grad[i][j][k][
l][j] * strain_tensor[q_point][k][
l]
3379 * dual_weights_cell_values[q_point](i) *
3380 * fe_values.JxW(q_point);
3390 * -------------------------------------------------------
3391 * compute face_integrals
3394 *
for (
unsigned int face_no=0;
3395 * face_no<GeometryInfo<dim>::faces_per_cell;
3398 *
if (cell->face(face_no)->at_boundary())
3400 *
for (
unsigned int id=0;
id!=dim; ++id)
3402 * face_integrals[cell->face(face_no)](id) = 0;
3407 *
if ((cell->neighbor(face_no)->has_children() ==
false) &&
3408 * (cell->neighbor(face_no)->level() == cell->level()) &&
3409 * (cell->neighbor(face_no)->index() < cell->index()))
3412 *
if (cell->at_boundary(face_no) ==
false)
3413 *
if (cell->neighbor(face_no)->level() < cell->level())
3417 *
if (cell->face(face_no)->has_children() ==
false)
3421 * ------------- integrate_over_regular_face -----------
3424 * fe_face_values_cell.reinit(cell, face_no);
3425 * fe_face_values_cell.get_function_gradients (primal_solution,
3429 * ExcInternalError());
3430 *
const unsigned int
3431 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3433 * neighbor = cell->neighbor(face_no);
3435 * fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
3436 * fe_face_values_neighbor.get_function_gradients (primal_solution,
3439 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3441 * q_cell_strain_tensor = 0.;
3442 * q_neighbor_strain_tensor = 0.;
3443 *
for (
unsigned int i=0; i!=dim; ++i)
3445 *
for (
unsigned int j=0; j!=dim; ++j)
3447 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3448 * cell_grads[q_point][j][i] );
3449 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3450 * neighbor_grads[q_point][j][i] );
3454 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3455 * cell_stress_strain_tensor);
3456 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3457 * neighbor_stress_strain_tensor);
3459 * jump_residual[q_point] = 0.;
3460 *
for (
unsigned int i=0; i!=dim; ++i)
3462 *
for (
unsigned int j=0; j!=dim; ++j)
3464 *
for (
unsigned int k=0; k!=dim; ++k)
3466 *
for (
unsigned int l=0;
l!=dim; ++
l)
3468 * jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
3469 * q_cell_strain_tensor[k][l]
3471 * neighbor_stress_strain_tensor[i][j][k][l]*
3472 * q_neighbor_strain_tensor[k][l] )*
3473 * fe_face_values_cell.normal_vector(q_point)[j];
3481 * fe_face_values_cell.get_function_values (dual_weights,
3482 * dual_weights_face_values);
3485 * face_integral_vector = 0;
3486 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3488 *
for (
unsigned int i=0; i!=dim; ++i)
3490 * face_integral_vector(i) += jump_residual[q_point](i) *
3491 * dual_weights_face_values[q_point](i) *
3492 * fe_face_values_cell.JxW(q_point);
3496 *
Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
3497 * ExcInternalError());
3499 *
for (
unsigned int i=0; i!=dim; ++i)
3501 *
Assert (face_integrals[cell->face(face_no)](i) == -1e20,
3502 * ExcInternalError());
3503 * face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
3509 * -----------------------------------------------------
3517 * ------------- integrate_over_irregular_face ---------
3521 * face = cell->face(face_no);
3523 * neighbor = cell->neighbor(face_no);
3525 * ExcInternalError());
3526 *
Assert (neighbor->has_children(),
3527 * ExcInternalError());
3529 *
const unsigned int
3530 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3532 *
for (
unsigned int subface_no=0;
3533 * subface_no<face->n_children(); ++subface_no)
3536 * neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
3537 *
Assert (neighbor_child->face(neighbor_neighbor) ==
3538 * cell->face(face_no)->child(subface_no),
3539 * ExcInternalError());
3541 * fe_subface_values_cell.reinit (cell, face_no, subface_no);
3542 * fe_subface_values_cell.get_function_gradients (primal_solution,
3544 * fe_face_values_neighbor.reinit (neighbor_child,
3545 * neighbor_neighbor);
3546 * fe_face_values_neighbor.get_function_gradients (primal_solution,
3549 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3551 * q_cell_strain_tensor = 0.;
3552 * q_neighbor_strain_tensor = 0.;
3553 *
for (
unsigned int i=0; i!=dim; ++i)
3555 *
for (
unsigned int j=0; j!=dim; ++j)
3557 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3558 * cell_grads[q_point][j][i] );
3559 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3560 * neighbor_grads[q_point][j][i] );
3564 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3565 * cell_stress_strain_tensor);
3566 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3567 * neighbor_stress_strain_tensor);
3569 * jump_residual[q_point] = 0.;
3570 *
for (
unsigned int i=0; i!=dim; ++i)
3572 *
for (
unsigned int j=0; j!=dim; ++j)
3574 *
for (
unsigned int k=0; k!=dim; ++k)
3576 *
for (
unsigned int l=0;
l!=dim; ++
l)
3578 * jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
3579 * q_cell_strain_tensor[k][l]
3581 * neighbor_stress_strain_tensor[i][j][k][l]*
3582 * q_neighbor_strain_tensor[k][l] )*
3583 * fe_face_values_neighbor.normal_vector(q_point)[j];
3591 * fe_face_values_neighbor.get_function_values (dual_weights,
3592 * dual_weights_face_values);
3595 * face_integral_vector = 0;
3596 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3598 *
for (
unsigned int i=0; i!=dim; ++i)
3600 * face_integral_vector(i) += jump_residual[q_point](i) *
3601 * dual_weights_face_values[q_point](i) *
3602 * fe_face_values_neighbor.JxW(q_point);
3606 *
for (
unsigned int i=0; i!=dim; ++i)
3608 * face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
3615 *
for (
unsigned int subface_no=0;
3616 * subface_no<face->n_children(); ++subface_no)
3618 *
Assert (face_integrals.find(face->child(subface_no)) !=
3619 * face_integrals.end(),
3620 * ExcInternalError());
3621 *
for (
unsigned int i=0; i!=dim; ++i)
3623 *
Assert (face_integrals[face->child(subface_no)](i) != -1e20,
3624 * ExcInternalError());
3625 *
sum(i) += face_integrals[face->child(subface_no)](i);
3628 *
for (
unsigned int i=0; i!=dim; ++i)
3630 * face_integrals[face](i) =
sum(i);
3636 * -----------------------------------------------------
3646 * ----------------------------------------------------------
3653 * cell = dof_handler_dual.begin_active();
3654 *
for (; cell!=endc; ++cell, ++present_cell)
3655 *
if (cell->is_locally_owned())
3657 *
for (
unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
3660 *
Assert(face_integrals.find(cell->face(face_no)) !=
3661 * face_integrals.end(),
3662 * ExcInternalError());
3664 *
for (
unsigned int id=0;
id!=dim; ++id)
3666 * error_indicators_vector[present_cell](id)
3667 * -= 0.5*face_integrals[cell->face(face_no)](id);
3672 * estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
3682 * <a name=
"ThecodePlasticityContactProblemcodeclasstemplate"></a>
3683 * <h3>The <code>PlasticityContactProblem</code>
class template</h3>
3687 * This is the main
class of this program and supplies all
functions
3688 * and variables needed to describe
3689 * the nonlinear contact problem. It is
3690 * close to @ref step_41
"step-41" but with some additional
3691 * features like handling hanging nodes,
3692 * a Newton method,
using Trilinos and p4est
3693 *
for parallel distributed computing.
3694 * To deal with hanging nodes makes
3695 * life a bit more complicated since
3697 * We create a Newton method
for the
3698 * active
set method
for the contact
3699 * situation and to handle the nonlinear
3700 *
operator for the constitutive law.
3704 * The
general layout of
this class is very much like for most other tutorial programs.
3705 * To make our life a bit easier,
this class reads a
set of input parameters from an input file. These
3706 * parameters,
using the
ParameterHandler class, are declared in the <code>declare_parameters</code>
3707 * function (which is
static so that it can be called before we even create an
object of the current
3708 * type), and a
ParameterHandler object that has been used to read an input file will then be passed
3709 * to the constructor of
this class.
3713 * The remaining member
functions are by and large as we have seen in several of the other tutorial
3714 * programs, though with additions
for the current nonlinear system. We will comment on their purpose
3715 * as we get to them further below.
3718 *
template <
int dim>
3719 *
class ElastoPlasticProblem
3729 *
void make_grid ();
3730 *
void setup_system ();
3731 *
void compute_dirichlet_constraints ();
3735 *
void solve_newton_system ();
3736 *
void solve_newton ();
3737 *
void compute_error ();
3739 *
void refine_grid ();
3741 *
void output_results (
const std::string &filename_base);
3745 * Next are three
functions that handle the history variables stored in each
3746 * quadrature
point. The
first one is called before the
first timestep to
3747 *
set up a pristine state
for the history variables. It only works on
3748 * those quadrature points on cells that belong to the present processor:
3751 *
void setup_quadrature_point_history ();
3755 * The
second one updates the history variables at the
end of each
3759 *
void update_quadrature_point_history ();
3763 * As far as member variables are concerned, we start with ones that we use to
3764 * indicate the MPI universe
this program runs on, and then two
numbers
3765 * telling us how many participating processors there are, and where in
3766 *
this world we are., a stream we use to let
3767 * exactly one processor produce output to the console (see @ref step_17
"step-17") and
3768 * a variable that is used to time the various sections of the program:
3779 * The next
group describes the mesh and the finite element space.
3780 * In particular,
for this parallel program, the finite element
3781 * space has associated with it variables that indicate which degrees
3782 * of freedom live on the current processor (the index sets, see
3783 * also @ref step_40
"step-40" and the @ref distributed documentation module) as
3784 * well as a variety of constraints: those imposed by hanging nodes,
3785 * by Dirichlet boundary conditions, and by the active
set of
3787 * here, the
first only contains hanging node constraints, the
3788 *
second also those associated with Dirichlet boundary conditions,
3789 * and the third these plus the contact constraints.
3793 * The variable <code>active_set</code> consists of those degrees
3794 * of freedom constrained by the contact, and we use
3795 * <code>fraction_of_plastic_q_points_per_cell</code> to keep
3796 * track of the fraction of quadrature points on each cell where
3797 * the stress equals the yield stress. The latter is only used to
3798 * create graphical output showing the plastic zone, but not
for
3799 * any further computation; the variable is a member variable of
3800 *
this class since the information is computed as a by-product
3801 * of computing the residual, but is used only much later. (Note
3802 * that the vector is a vector of length
equal to the number of
3803 * active cells on the <i>local mesh</i>; it is never used to
3804 * exchange information between processors and can therefore be
3805 * a regular deal.II vector.)
3808 *
const unsigned int n_initial_global_refinements;
3811 *
const unsigned int fe_degree;
3825 * One difference of
this program is that we declare the quadrature
3826 * formula in the
class declaration. The reason is that in all the other
3827 * programs, it didn
't do much harm if we had used different quadrature
3828 * formulas when computing the matrix and the right hand side, for
3829 * example. However, in the present case it does: we store information in
3830 * the quadrature points, so we have to make sure all parts of the program
3831 * agree on where they are and how many there are on each cell. Thus, let
3832 * us first declare the quadrature formula that will be used throughout...
3835 * const QGauss<dim> quadrature_formula;
3836 * const QGauss<dim - 1> face_quadrature_formula;
3840 * ... and then also have a vector of history objects, one per quadrature
3841 * point on those cells for which we are responsible (i.e. we don't store
3842 * history data
for quadrature points on cells that are owned by other
3846 * std::vector<PointHistory<dim> > quadrature_point_history;
3850 * The way
this object is accessed is through a <code>user pointer</code>
3851 * that each cell, face, or edge holds: it is a <code>
void*</code> pointer
3852 * that can be used by application programs to associate arbitrary data to
3853 * cells, faces, or edges. What the program actually does with
this data
3854 * is within its own responsibility, the library just allocates some space
3855 *
for these pointers, and application programs can
set and read the
3856 * pointers
for each of these objects.
3863 * The next block of variables corresponds to the solution
3864 * and the linear systems we need to form. In particular,
this
3865 * includes the Newton
matrix and right hand side; the vector
3866 * that corresponds to the residual (i.e., the Newton right hand
3867 * side) but from which we have not eliminated the various
3868 * constraints and that is used to determine which degrees of
3869 * freedom need to be constrained in the next iteration; and
3870 * a vector that corresponds to the
diagonal of the @f$B@f$
matrix
3871 * briefly mentioned in the introduction and discussed in the
3872 * accompanying paper.
3884 * The next block of variables is then related to the time dependent
3885 * nature of the problem: they denote the length of the time interval
3886 * which we want to simulate, the present time and number of time step,
3887 * and length of present timestep:
3890 *
double present_time;
3891 *
double present_timestep;
3893 *
unsigned int timestep_no;
3897 * The next block contains the variables that describe the material
3901 *
const double e_modulus, nu, sigma_0,
gamma;
3902 * ConstitutiveLaw<dim> constitutive_law;
3906 * And then there is an assortment of other variables that are used
3907 * to identify the mesh we are asked to build as selected by the
3908 * parameter file, the obstacle that is being pushed into the
3909 * deformable body, the mesh refinement strategy, whether to transfer
3910 * the solution from one mesh to the next, and how many mesh
3911 * refinement cycles to perform. As possible, we mark these kinds
3912 * of variables as <code>
const</code> to help the reader identify
3913 * which ones may or may not be modified later on (the output directory
3914 * being an exception -- it is never modified outside the constructor
3915 * but it is awkward to initialize in the member-initializer-list
3916 * following the colon in the constructor since there we have only
3917 * one shot at setting it; the same is
true for the mesh refinement
3921 * const
std::string base_mesh;
3923 *
struct RefinementStrategy
3928 * refine_percentage,
3932 *
typename RefinementStrategy::value refinement_strategy;
3934 *
struct ErrorEstimationStrategy
3940 * weighted_residual_error,
3941 * weighted_kelly_error
3944 *
typename ErrorEstimationStrategy::value error_estimation_strategy;
3948 *
const bool transfer_solution;
3949 * std::string output_dir;
3954 *
unsigned int current_refinement_cycle;
3956 *
const double max_relative_error;
3957 *
float relative_error;
3959 *
const bool show_stresses;
3966 * <a name=
"ImplementationofthecodePlasticityContactProblemcodeclass"></a>
3967 * <h3>Implementation of the <code>PlasticityContactProblem</code>
class</h3>
3972 * <a name=
"PlasticityContactProblemdeclare_parameters"></a>
3973 * <h4>PlasticityContactProblem::declare_parameters</h4>
3977 * Let us start with the declaration of
run-time parameters that can be
3978 * selected in the input file. These
values will be read back in the
3979 * constructor of
this class to initialize the member variables of this
3983 *
template <
int dim>
3987 * prm.declare_entry(
"polynomial degree",
"1",
3989 *
"Polynomial degree of the FE_Q finite element space, typically 1 or 2.");
3990 * prm.declare_entry(
"number of initial refinements",
"2",
3992 *
"Number of initial global mesh refinement steps before "
3993 *
"the first computation.");
3994 * prm.declare_entry(
"refinement strategy",
"percentage",
3996 *
"Mesh refinement strategy:\n"
3997 *
" global: one global refinement\n"
3998 *
" percentage: a fixed percentage of cells gets refined using the selected error estimator.");
3999 * prm.declare_entry(
"error estimation strategy",
"kelly_error",
4001 *
"Error estimation strategy:\n"
4002 *
" kelly_error: Kelly error estimator\n"
4003 *
" residual_error: residual-based error estimator\n"
4004 *
" weighted_residual_error: dual weighted residual (Goal-oriented) error estimator.\n");
4005 * prm.declare_entry(
"maximum relative error",
"0.05",
4007 *
"maximum relative error which plays the role of a criteria for refinement.");
4008 * prm.declare_entry(
"number of cycles",
"5",
4010 *
"Number of adaptive mesh refinement cycles to run.");
4011 * prm.declare_entry(
"output directory",
"",
4013 *
"Directory for output files (graphical output and benchmark "
4014 *
"statistics). If empty, use the current directory.");
4015 * prm.declare_entry(
"transfer solution",
"true",
4017 *
"Whether the solution should be used as a starting guess "
4018 *
"for the next finer mesh. If false, then the iteration starts at "
4019 *
"zero on every mesh.");
4020 * prm.declare_entry(
"base mesh",
"Thick_tube_internal_pressure",
4022 *
"Perforated_strip_tension|Cantiliver_beam_3d"),
4023 *
"Select the shape of the domain: 'box' or 'half sphere'");
4024 * prm.declare_entry(
"elasticity modulus",
"2.e11",
4026 *
"Elasticity modulus of the material in MPa (N/mm2)");
4027 * prm.declare_entry(
"Poissons ratio",
"0.3",
4029 *
"Poisson's ratio of the material");
4030 * prm.declare_entry(
"yield stress",
"2.e11",
4032 *
"Yield stress of the material in MPa (N/mm2)");
4033 * prm.declare_entry(
"isotropic hardening parameter",
"0.",
4035 *
"Isotropic hardening parameter of the material");
4036 * prm.declare_entry(
"show stresses",
"false",
4038 *
"Whether illustrates the stresses and von Mises stresses or not.");
4047 * <a name=
"ThecodePlasticityContactProblemcodeconstructor"></a>
4048 * <h4>The <code>PlasticityContactProblem</code> constructor</h4>
4052 * Given the declarations of member variables as well as the
4053 * declarations of
run-time parameters that are read from the input
4054 * file, there is
nothing surprising in
this constructor. In the body
4055 * we initialize the mesh refinement strategy and the output directory,
4056 * creating such a directory
if necessary.
4059 *
template <
int dim>
4060 * ElastoPlasticProblem<dim>::
4063 * mpi_communicator(MPI_COMM_WORLD),
4066 * pcout(std::cout, this_mpi_process == 0),
4070 * n_initial_global_refinements (prm.get_integer(
"number of initial refinements")),
4072 * fe_degree (prm.get_integer(
"polynomial degree")),
4075 * quadrature_formula (fe_degree + 1),
4076 * face_quadrature_formula (fe_degree + 1),
4078 * e_modulus (prm.get_double(
"elasticity modulus")),
4079 * nu (prm.get_double(
"Poissons ratio")),
4080 * sigma_0(prm.get_double(
"yield stress")),
4081 *
gamma (prm.get_double(
"isotropic hardening parameter")),
4082 * constitutive_law (e_modulus,
4087 * base_mesh (prm.get(
"base mesh")),
4089 * transfer_solution (prm.get_bool(
"transfer solution")),
4091 * table_results_2(),
4092 * table_results_3(),
4093 * max_relative_error (prm.get_double(
"maximum relative error")),
4094 * show_stresses (prm.get_bool(
"show stresses"))
4096 * std::string strat = prm.get(
"refinement strategy");
4097 *
if (strat ==
"global")
4098 * refinement_strategy = RefinementStrategy::refine_global;
4099 *
else if (strat ==
"percentage")
4100 * refinement_strategy = RefinementStrategy::refine_percentage;
4104 * strat = prm.get(
"error estimation strategy");
4105 *
if (strat ==
"kelly_error")
4106 * error_estimation_strategy = ErrorEstimationStrategy::kelly_error;
4107 *
else if (strat ==
"residual_error")
4108 * error_estimation_strategy = ErrorEstimationStrategy::residual_error;
4109 *
else if (strat ==
"weighted_residual_error")
4110 * error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error;
4114 * output_dir = prm.get(
"output directory");
4115 *
if (output_dir !=
"" && *(output_dir.rbegin()) !=
'/')
4116 * output_dir +=
"/";
4117 * mkdir(output_dir.c_str(), 0777);
4119 * pcout <<
" Using output directory '" << output_dir <<
"'" << std::endl;
4120 * pcout <<
" FE degree " << fe_degree << std::endl;
4121 * pcout <<
" transfer solution "
4122 * << (transfer_solution ?
"true" :
"false") << std::endl;
4130 * <a name=
"PlasticityContactProblemmake_grid"></a>
4131 * <h4>PlasticityContactProblem::make_grid</h4>
4135 * The next block deals with constructing the starting mesh.
4136 * We will use the following helper function and the
first
4137 * block of the <code>make_grid()</code> to construct a
4138 * mesh that corresponds to a half sphere. deal.II has a function
4139 * that creates such a mesh, but it is in the wrong location
4140 * and facing the wrong direction, so we need to
shift and
rotate
4141 * it a bit before
using it.
4145 * For later reference, as described in the documentation of
4147 * has boundary indicator zero,
while the remainder has boundary
4152 * rotate_half_sphere (
const Point<3> &in)
4154 *
return Point<3>(in(2), in(1), -in(0));
4157 *
template <
int dim>
4159 * ElastoPlasticProblem<dim>::make_grid ()
4161 *
if (base_mesh ==
"Timoshenko beam")
4165 *
const double length = .48,
4169 * point_2(length,
depth/2);
4171 * std::vector<unsigned int> repetitions(2);
4172 * repetitions[0] = 4;
4173 * repetitions[1] = 1;
4179 * give the indicators to boundaries
for specification,
4186 * |________________|
4188 * 0 to essential boundary conditions (left edge) which are as
default
4189 * 100 to the null boundaries (upper and lower edges) where we
do not need to take care of them
4190 * 5 to the natural boundaries (right edge)
for imposing the traction force
4196 *
for (; cell!=endc; ++cell)
4200 *
if ( std::fabs(cell->face(face)->center()(0)-length) < 1e-12 )
4202 * cell->face(face)->set_manifold_id(5);
4204 *
else if ( ( std::fabs(cell->face(face)->center()(1)-(
depth/2)) < 1e-12 )
4206 * ( std::fabs(cell->face(face)->center()(1)-(-
depth/2)) < 1e-12 ) )
4208 * cell->face(face)->set_manifold_id(100);
4214 *
triangulation.refine_global(n_initial_global_refinements);
4217 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4221 * Example 1 from the paper: Zhong Z., .... A
new numerical method
for determining
4222 * collapse load-carrying capacity of structure made of elasto-plastic material,
4223 * J. Cent. South Univ. (2014) 21: 398-404
4229 *
const double inner_radius = .1,
4230 * outer_radius = .2;
4232 *
center, inner_radius, outer_radius,
4237 * give the indicators to boundaries
for specification,
4255 * 0 - inner boundary - natural boundary condition - impose the traction force
4256 * 1 - outer boundary -
free boundary - we
do not need to take care of them
4257 * 2 - left boundary - essential boundary condition - constrained to move along the x direction
4258 * 3 - bottom boundary - essential boundary condition - constrained to move along the y direction
4265 *
triangulation.set_manifold (0, inner_boundary_description);
4268 *
triangulation.set_manifold (1, outer_boundary_description);
4270 *
triangulation.refine_global(n_initial_global_refinements);
4276 *
else if (base_mesh ==
"Perforated_strip_tension")
4280 * Example 2 from the paper: Zhong Z., .... A
new numerical method
for determining
4281 * collapse load-carrying capacity of structure made of elasto-plastic material,
4282 * J. Cent. South Univ. (2014) 21: 398-404
4287 *
const int dim_2d = 2;
4289 *
const double inner_radius = 0.05,
4290 * outer_radius = 0.1,
4292 * thickness = 0.004;
4305 *
const double eps = 1
e-7 * inner_radius;
4310 * center_2d, inner_radius, outer_radius,
4315 * Modify the triangulation_1
4319 * cell = triangulation_1.begin_active(),
4320 * endc = triangulation_1.end();
4321 * std::vector<bool> treated_vertices(triangulation_1.n_vertices(),
false);
4322 *
for (; cell != endc; ++cell)
4324 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4325 *
if (cell->face(f)->at_boundary() && cell->face(f)->center()(0)>eps &&
4326 * cell->face(f)->center()(1)>eps )
4333 *
point(0) = cell->face(f)->center()(0) - center_2d(0);
4334 *
point(1) = cell->face(f)->center()(1) - center_2d(1);
4337 *
for (
unsigned int v=0; v < GeometryInfo<dim_2d>::vertices_per_face; ++v)
4339 *
unsigned int vv = cell->face(f)->vertex_index(v);
4340 *
if (treated_vertices[vv] ==
false)
4342 * treated_vertices[vv] =
true;
4345 * cell->face(f)->vertex(v) = center_2d+
Point<dim_2d>(outer_radius,outer_radius);
4358 * Make the triangulation_2, a rectangular above the triangulation_1
4363 * point2 (outer_radius, height);
4371 * make the triangulation_2d and
refine it
4377 * Merge the two triangulation_1 and triangulation_2
4384 * Assign boundary indicators to the boundary faces
4405 * cell = triangulation_2d.begin_active(),
4406 * endc = triangulation_2d.end();
4407 *
for (; cell != endc; ++cell)
4409 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4411 *
if (cell->face(f)->at_boundary())
4413 *
if ( std::fabs(cell->face(f)->center()(1)) < eps )
4415 * cell->face(f)->set_manifold_id(1);
4417 *
else if ( std::fabs(cell->face(f)->center()(0)-outer_radius) < eps )
4419 * cell->face(f)->set_manifold_id(2);
4421 *
else if ( std::fabs(cell->face(f)->center()(1)-height) < eps )
4423 * cell->face(f)->set_manifold_id(3);
4425 *
else if ( std::fabs(cell->face(f)->center()(0)) < eps )
4427 * cell->face(f)->set_manifold_id(4);
4431 * cell->face(f)->set_all_boundary_ids(10);
4441 * triangulation_2d.set_manifold (10, inner_boundary_description);
4443 * triangulation_2d.refine_global(3);
4445 * triangulation_2d.reset_manifold (10);
4450 * Extrude the triangulation_2d and make it 3
d
4460 * Assign boundary indicators to the boundary faces
4486 *
for (; cell != endc; ++cell)
4488 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4490 *
if (cell->face(f)->at_boundary())
4492 * dist_vector = cell->face(f)->center() -
center;
4494 *
if ( std::fabs(dist_vector[1]) < eps )
4496 * cell->face(f)->set_manifold_id(1);
4498 *
else if ( std::fabs(dist_vector[0]-outer_radius) < eps )
4500 * cell->face(f)->set_manifold_id(2);
4502 *
else if ( std::fabs(dist_vector[1]-height) < eps )
4504 * cell->face(f)->set_manifold_id(3);
4506 *
else if ( std::fabs(dist_vector[0]) < eps )
4508 * cell->face(f)->set_manifold_id(4);
4510 *
else if ( std::fabs(dist_vector[2]) < eps )
4512 * cell->face(f)->set_manifold_id(5);
4514 *
else if ( std::fabs(dist_vector[2]-thickness) < eps )
4516 * cell->face(f)->set_manifold_id(6);
4520 * cell->face(f)->set_all_boundary_ids(10);
4530 *
triangulation.set_manifold (10, inner_boundary_description);
4532 *
triangulation.refine_global(n_initial_global_refinements);
4537 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4541 * A rectangular tube made of Aluminium
4543 * approximation of beam 17250
4544 * units are in meter
4552 *
const int dim_2d = 2;
4554 *
const double length = .7,
4557 * thickness_web = 10
e-3,
4558 * thickness_flange = 10
e-3;
4566 *
const double eps = 1
e-7 * width;
4569 * Make the triangulation_b, a rectangular at the bottom of rectangular tube
4574 * point2 (width/2, -(height/2)+thickness_flange);
4576 * std::vector<unsigned int> repetitions(dim_2d);
4577 * repetitions[0] = 8;
4578 * repetitions[1] = 1;
4585 * Make the triangulation_t, a rectangular at the top of rectangular tube
4589 *
const Point<dim_2d> point1 (-width/2, (height/2)-thickness_flange),
4590 * point2 (width/2, height/2);
4592 * std::vector<unsigned int> repetitions(dim_2d);
4593 * repetitions[0] = 8;
4594 * repetitions[1] = 1;
4601 * Make the triangulation_l, a rectangular at the left of rectangular tube
4605 *
const Point<dim_2d> point1 (-width/2, -(height/2)+thickness_flange),
4606 * point2 (-(width/2)+thickness_web, (height/2)-thickness_flange);
4608 * std::vector<unsigned int> repetitions(dim_2d);
4609 * repetitions[0] = 1;
4610 * repetitions[1] = 18;
4617 * Make the triangulation_r, a rectangular at the right of rectangular tube
4621 *
const Point<dim_2d> point1 ((width/2)-thickness_web, -(height/2)+thickness_flange),
4622 * point2 (width/2, (height/2)-thickness_flange);
4624 * std::vector<unsigned int> repetitions(dim_2d);
4625 * repetitions[0] = 1;
4626 * repetitions[1] = 18;
4633 * make the triangulation_2d
4639 * merging every two triangles to make triangulation_2d
4643 * triangulation_blr;
4652 * Extrude the triangulation_2d and make it 3
d
4655 *
const unsigned int n_slices =
static_cast<int>(length*1000/20) + 1;
4661 * Assign boundary indicators to the boundary faces
4705 *
for (; cell != endc; ++cell)
4707 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4709 *
if (cell->face(f)->at_boundary())
4711 * dist_vector = cell->face(f)->center() -
center;
4713 *
if ( std::fabs(dist_vector[2]) < eps )
4715 * cell->face(f)->set_manifold_id(1);
4717 *
else if ( std::fabs(dist_vector[1]-(height/2)) < eps )
4719 * cell->face(f)->set_manifold_id(2);
4723 * cell->face(f)->set_all_boundary_ids(0);
4732 *
triangulation.refine_global(n_initial_global_refinements);
4740 * pcout <<
" Number of active cells: "
4750 * <a name=
"PlasticityContactProblemsetup_system"></a>
4751 * <h4>PlasticityContactProblem::setup_system</h4>
4755 * The next piece in the puzzle is to
set up the
DoFHandler, resize
4756 * vectors and take care of various other status variables such as
4757 *
index sets and constraint matrices.
4761 * In the following, each
group of operations is put into a brace-enclosed
4762 * block that is being timed by the variable declared at the top of the
4764 * timed section, the destructor that is called at the end of the block
4768 *
template <
int dim>
4770 * ElastoPlasticProblem<dim>::setup_system ()
4776 * dof_handler.distribute_dofs(fe);
4777 * pcout <<
" Number of degrees of freedom: "
4778 * << dof_handler.n_dofs()
4781 * locally_owned_dofs = dof_handler.locally_owned_dofs();
4782 * locally_relevant_dofs.clear();
4784 * locally_relevant_dofs);
4790 * constraints_hanging_nodes.reinit(locally_relevant_dofs);
4792 * constraints_hanging_nodes);
4793 * constraints_hanging_nodes.close();
4795 * pcout <<
" Number of active cells: "
4797 * <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
4800 * compute_dirichlet_constraints();
4806 *
if (timestep_no==1 || current_refinement_cycle!=0)
4808 * solution.reinit(locally_relevant_dofs, mpi_communicator);
4810 * incremental_displacement.reinit(locally_relevant_dofs, mpi_communicator);
4811 * newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
4812 * newton_rhs_residual.reinit(locally_owned_dofs, mpi_communicator);
4813 * fraction_of_plastic_q_points_per_cell.reinit(
triangulation.n_active_cells());
4818 * Finally, we
set up sparsity patterns and matrices.
4819 * We temporarily (ab)use the system
matrix to also build the (diagonal)
4820 *
matrix that we use in eliminating degrees of freedom that are in contact
4821 * with the obstacle, but we then immediately
set the Newton
matrix back
4828 * mpi_communicator);
4831 * constraints_dirichlet_and_hanging_nodes,
false,
4832 * this_mpi_process);
4834 * newton_matrix.reinit(sp);
4842 * <a name=
"PlasticityContactProblemcompute_dirichlet_constraints"></a>
4843 * <h4>PlasticityContactProblem::compute_dirichlet_constraints</h4>
4847 * This function, broken out of the preceding one, computes the constraints
4848 * associated with Dirichlet-type boundary conditions and puts them into the
4849 * <code>constraints_dirichlet_and_hanging_nodes</code> variable by merging
4850 * with the constraints that come from hanging nodes.
4854 * As laid out in the introduction, we need to distinguish between two
4856 * - If the domain is a box, we
set the displacement to zero at the bottom,
4857 * and allow vertical movement in z-direction along the sides. As
4858 * shown in the <code>make_grid()</code> function, the former corresponds
4859 * to boundary indicator 6, the latter to 8.
4860 * - If the domain is a half sphere, then we impose zero displacement along
4861 * the curved part of the boundary, associated with boundary indicator zero.
4864 *
template <
int dim>
4866 * ElastoPlasticProblem<dim>::compute_dirichlet_constraints ()
4868 * constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
4869 * constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
4871 * std::vector<bool> component_mask(dim);
4873 *
if (base_mesh ==
"Timoshenko beam")
4877 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4878 * constraints_dirichlet_and_hanging_nodes,
4881 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4885 * the boundary x = 0
4888 * component_mask[0] =
true;
4889 * component_mask[1] =
false;
4892 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4893 * constraints_dirichlet_and_hanging_nodes,
4897 * the boundary y = 0
4900 * component_mask[0] =
false;
4901 * component_mask[1] =
true;
4904 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4905 * constraints_dirichlet_and_hanging_nodes,
4908 *
else if (base_mesh ==
"Perforated_strip_tension")
4912 * the boundary x = 0
4915 * component_mask[0] =
true;
4916 * component_mask[1] =
false;
4917 * component_mask[2] =
false;
4920 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4921 * constraints_dirichlet_and_hanging_nodes,
4925 * the boundary y = 0
4928 * component_mask[0] =
false;
4929 * component_mask[1] =
true;
4930 * component_mask[2] =
false;
4933 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4934 * constraints_dirichlet_and_hanging_nodes,
4938 * the boundary y = imposed incremental displacement
4941 * component_mask[0] =
false;
4942 * component_mask[1] =
true;
4943 * component_mask[2] =
false;
4946 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4947 * constraints_dirichlet_and_hanging_nodes,
4950 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4954 * the boundary x = y = z = 0
4957 * component_mask[0] =
true;
4958 * component_mask[1] =
true;
4959 * component_mask[2] =
true;
4962 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4963 * constraints_dirichlet_and_hanging_nodes,
4972 * constraints_dirichlet_and_hanging_nodes.close();
4979 * <a name=
"PlasticityContactProblemassemble_newton_system"></a>
4980 * <h4>PlasticityContactProblem::assemble_newton_system</h4>
4984 * Given the complexity of the problem, it may come as a bit of a surprise
4985 * that assembling the linear system we have to solve in each Newton iteration
4986 * is actually fairly straightforward. The following function builds the Newton
4987 * right hand side and Newton
matrix. It looks fairly innocent because the
4988 * heavy lifting happens in the
call to
4989 * <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
4991 * constraints we have previously computed.
4994 *
template <
int dim>
4996 * ElastoPlasticProblem<dim>::
5003 *
if (base_mesh ==
"Timoshenko beam")
5005 * traction_surface_id = 5;
5007 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5009 * traction_surface_id = 0;
5011 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5013 * traction_surface_id = 2;
5027 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5028 *
const unsigned int n_q_points = quadrature_formula.size();
5029 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5032 *
const EquationData::BodyForce<dim> body_force;
5033 * std::vector<Vector<double> > body_force_values(n_q_points,
5036 *
const EquationData::
5037 * IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5038 * std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5044 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5048 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5051 * std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
5054 * cell = dof_handler.begin_active(),
5055 * endc = dof_handler.end();
5059 *
for (; cell != endc; ++cell)
5060 *
if (cell->is_locally_owned())
5062 * fe_values.reinit(cell);
5066 * fe_values[displacement].get_function_symmetric_gradients(delta_linearization_point,
5067 * incremental_strain_tensor);
5071 * For assembling the local right hand side contributions, we need
5072 * to access the prior linearized stress
value in
this quadrature
5073 *
point. To get it, we use the user pointer of
this cell that
5074 * points into the global array to the quadrature
point data
5075 * corresponding to the
first quadrature
point of the present cell,
5076 * and then add an offset corresponding to the
index of the
5077 * quadrature
point we presently consider:
5080 *
const PointHistory<dim> *local_quadrature_points_history
5081 * =
reinterpret_cast<PointHistory<dim>*
>(cell->user_pointer());
5082 *
Assert (local_quadrature_points_history >=
5083 * &quadrature_point_history.front(),
5084 * ExcInternalError());
5085 *
Assert (local_quadrature_points_history <
5086 * &quadrature_point_history.back(),
5087 * ExcInternalError());
5091 * In addition, we need the
values of the external body forces at
5092 * the quadrature points on
this cell:
5095 * body_force.vector_value_list(fe_values.get_quadrature_points(),
5096 * body_force_values);
5098 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5101 * tmp_strain_tensor_qpoint = local_quadrature_points_history[q_point].old_strain
5102 * + incremental_strain_tensor[q_point];
5106 * constitutive_law.get_linearized_stress_strain_tensors(tmp_strain_tensor_qpoint,
5107 * stress_strain_tensor_linearized,
5108 * stress_strain_tensor);
5111 *
for (
unsigned int i = 0; i < dim; ++i)
5113 * rhs_values_body_force[i] = body_force_values[q_point][i];
5116 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5120 * Having computed the stress-strain tensor and its linearization,
5121 * we can now put together the parts of the
matrix and right hand side.
5122 * In both, we need the linearized stress-strain tensor times the
5123 *
symmetric gradient of @f$\varphi_i@f$, i.e. the term @f$I_\Pi\varepsilon(\varphi_i)@f$,
5124 * so we introduce an abbreviation of
this term. Recall that the
5125 *
matrix corresponds to the bilinear form
5126 * @f$A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))@f$ in the
5127 * notation of the accompanying publication, whereas the right
5128 * hand side is @f$F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))@f$
5129 * where @f$u@f$ is the current linearization points (typically the last solution).
5130 * This might suggest that the right hand side will be zero
if the material
5131 * is completely elastic (where @f$I_\Pi=P_\Pi@f$) but
this ignores the fact
5132 * that the right hand side will also contain contributions from
5133 * non-homogeneous constraints due to the contact.
5137 * The code block that follows
this adds contributions that are due to
5138 * boundary forces, should there be any.
5142 * stress_phi_i = stress_strain_tensor_linearized
5143 * * fe_values[displacement].symmetric_gradient(i, q_point);
5145 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
5147 * * fe_values[displacement].symmetric_gradient(j, q_point)
5148 * * fe_values.JxW(q_point));
5152 * * incremental_strain_tensor[q_point] )
5154 * ( ( stress_strain_tensor
5155 * * fe_values[displacement].symmetric_gradient(i, q_point))
5156 * * tmp_strain_tensor_qpoint )
5158 * ( fe_values[displacement].
value(i, q_point)
5159 * * rhs_values_body_force )
5160 * ) * fe_values.JxW(q_point);
5165 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
5166 *
if (cell->face(face)->at_boundary()
5168 * cell->face(face)->boundary_id() == traction_surface_id)
5170 * fe_values_face.reinit(cell, face);
5172 * boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5173 * boundary_force_values);
5175 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
5178 *
for (
unsigned int i = 0; i < dim; ++i)
5180 * rhs_values[i] = boundary_force_values[q_point][i];
5182 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5183 * cell_rhs(i) += (fe_values_face[displacement].value(i, q_point)
5185 * * fe_values_face.JxW(q_point));
5189 * cell->get_dof_indices(local_dof_indices);
5190 * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_matrix, cell_rhs,
5191 * local_dof_indices,
5207 * <a name=
"PlasticityContactProblemcompute_nonlinear_residual"></a>
5208 * <h4>PlasticityContactProblem::compute_nonlinear_residual</h4>
5212 * The following function computes the nonlinear residual of the equation
5213 * given the current solution (or any other linearization point). This
5214 * is needed in the linear search algorithm where we need to
try various
5215 * linear combinations of previous and current (trial) solution to
5216 * compute the (real, globalized) solution of the current Newton step.
5220 * That said, in a slight abuse of the name of the function, it actually
5221 * does significantly more. For example, it also computes the vector
5222 * that corresponds to the Newton residual but without eliminating
5223 * constrained degrees of freedom. We need
this vector to compute contact
5224 * forces and, ultimately, to compute the next active
set. Likewise, by
5225 * keeping track of how many quadrature points we encounter on each cell
5226 * that show plastic yielding, we also compute the
5227 * <code>fraction_of_plastic_q_points_per_cell</code> vector that we
5228 * can later output to visualize the plastic zone. In both of these cases,
5229 * the results are not necessary as part of the line search, and so we may
5230 * be wasting a small amount of time computing them. At the same time,
this
5231 * information appears as a natural by-product of what we need to
do here
5232 * anyway, and we want to collect it once at the
end of each Newton
5233 * step, so we may as well
do it here.
5237 * The actual implementation of
this function should be rather obvious:
5240 *
template <
int dim>
5242 * ElastoPlasticProblem<dim>::
5246 *
if (base_mesh ==
"Timoshenko beam")
5248 * traction_surface_id = 5;
5250 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5252 * traction_surface_id = 0;
5254 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5256 * traction_surface_id = 2;
5271 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5272 *
const unsigned int n_q_points = quadrature_formula.size();
5273 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5275 *
const EquationData::BodyForce<dim> body_force;
5276 * std::vector<Vector<double> > body_force_values(n_q_points,
5279 *
const EquationData::
5280 * IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5281 * std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5286 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5290 * newton_rhs_residual = 0;
5292 * fraction_of_plastic_q_points_per_cell = 0;
5295 * cell = dof_handler.begin_active(),
5296 * endc = dof_handler.end();
5297 *
unsigned int cell_number = 0;
5298 *
for (; cell != endc; ++cell, ++cell_number)
5299 *
if (cell->is_locally_owned())
5301 * fe_values.reinit(cell);
5304 * std::vector<SymmetricTensor<2, dim> > strain_tensors(n_q_points);
5305 * fe_values[displacement].get_function_symmetric_gradients(linearization_point,
5308 * body_force.vector_value_list(fe_values.get_quadrature_points(),
5309 * body_force_values);
5311 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5314 *
const bool q_point_is_plastic
5315 * = constitutive_law.get_stress_strain_tensor(strain_tensors[q_point],
5316 * stress_strain_tensor);
5317 *
if (q_point_is_plastic)
5318 * ++fraction_of_plastic_q_points_per_cell(cell_number);
5321 *
for (
unsigned int i = 0; i < dim; ++i)
5323 * rhs_values_body_force[i] = body_force_values[q_point][i];
5326 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5328 * cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5329 * * rhs_values_body_force
5331 * strain_tensors[q_point]
5332 * * stress_strain_tensor
5333 * * fe_values[displacement].symmetric_gradient(i, q_point)
5335 * * fe_values.JxW(q_point);
5339 * cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5341 * * fe_values.JxW(q_point));
5345 *
for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5346 *
if (cell->face(face)->at_boundary()
5347 * && cell->face(face)->boundary_id() == traction_surface_id)
5349 * fe_values_face.reinit(cell, face);
5351 * boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5352 * boundary_force_values);
5354 *
for (
unsigned int q_point = 0; q_point < n_face_q_points;
5358 *
for (
unsigned int i = 0; i < dim; ++i)
5360 * rhs_values[i] = boundary_force_values[q_point][i];
5362 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5363 * cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values
5364 * * fe_values_face.JxW(q_point));
5368 * cell->get_dof_indices(local_dof_indices);
5369 * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_rhs,
5370 * local_dof_indices,
5371 * newton_rhs_residual);
5375 * fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
5387 * <a name=
"PlasticityContactProblemsolve_newton_system"></a>
5388 * <h4>PlasticityContactProblem::solve_newton_system</h4>
5392 * The last piece before we can discuss the actual Newton iteration
5393 * on a single mesh is the solver
for the linear systems. There are
5394 * a couple of complications that slightly obscure the code, but
5395 * mostly it is just setup then solve. Among the complications are:
5399 * - For the hanging nodes we have to
apply
5401 * This is necessary
if a hanging node with solution
value @f$x_0@f$
5402 * has one neighbor with
value @f$x_1@f$ which is in contact with the
5403 * obstacle and one neighbor @f$x_2@f$ which is not in contact. Because
5404 * the update
for the former will be prescribed, the hanging node constraint
5405 * will have an inhomogeneity and will look like @f$x_0 = x_1/2 + \text{gap}/2@f$.
5406 * So the corresponding entries in the
5407 * ride-hang-side are non-zero with a
5410 * - Like in @ref step_40
"step-40", we need to shuffle between vectors that
do and
do
5411 *
do not have ghost elements when solving or
using the solution.
5415 * The rest of the function is similar to @ref step_40
"step-40" and
5416 * @ref step_41
"step-41" except that we use a BiCGStab solver
5417 * instead of CG. This is due to the fact that
for very small hardening
5418 * parameters @f$\gamma@f$, the linear system becomes almost semidefinite though
5419 * still
symmetric. BiCGStab appears to have an easier time with such linear
5423 *
template <
int dim>
5425 * ElastoPlasticProblem<dim>::solve_newton_system ()
5430 * distributed_solution = incremental_displacement;
5432 * constraints_hanging_nodes.set_zero(distributed_solution);
5433 * constraints_hanging_nodes.set_zero(newton_rhs);
5437 * ------- Solver Bicgstab --- Preconditioner AMG -------------------
5444 * std::vector<std::vector<bool> > constant_modes;
5454 * additional_data.
w_cycle =
false;
5461 * preconditioner.initialize(newton_matrix, additional_data);
5476 *
const double relative_accuracy = 1
e-2;
5477 *
const double solver_tolerance = relative_accuracy
5478 * * newton_matrix.residual(tmp, distributed_solution,
5484 * solver_tolerance);
5486 * solver.solve(newton_matrix, distributed_solution,
5487 * newton_rhs, preconditioner);
5491 * pcout <<
" Error: " << solver_control.initial_value()
5492 * <<
" -> " << solver_control.last_value() <<
" in "
5493 * << solver_control.last_step() <<
" Bicgstab iterations."
5499 * ------- Solver CG --- Preconditioner SSOR -------------------
5507 * preconditioner.initialize(newton_matrix, additional_data);
5517 *
const double relative_accuracy = 1
e-8;
5520 *
const double relative_accuracy = 1
e-2;
5521 *
const double solver_tolerance = relative_accuracy
5522 * * newton_matrix.residual(tmp, distributed_solution,
5528 * solver_tolerance);
5532 * solver_tolerance);
5534 * solver.solve(newton_matrix, distributed_solution,
5535 * newton_rhs, preconditioner);
5537 * pcout <<
" Error: " << solver_control.initial_value()
5538 * <<
" -> " << solver_control.last_value() <<
" in "
5539 * << solver_control.last_step() <<
" CG iterations."
5544 * ........................................................
5550 * constraints_dirichlet_and_hanging_nodes.distribute(distributed_solution);
5552 * incremental_displacement = distributed_solution;
5559 * <a name=
"PlasticityContactProblemsolve_newton"></a>
5560 * <h4>PlasticityContactProblem::solve_newton</h4>
5564 * This is,
finally, the function that implements the damped Newton method
5565 * on the current mesh. There are two nested loops: the outer
loop for the Newton
5566 * iteration and the inner
loop for the line search which
5567 * will be used only
if necessary. To obtain a good and reasonable
5568 * starting
value we solve an elastic problem in very
first Newton step on each
5569 * mesh (or only on the
first mesh
if we transfer solutions between meshes). We
5570 *
do so by setting the yield stress to an unreasonably large
value in these
5571 * iterations and then setting it back to the correct
value in subsequent
5576 * Other than
this, the top part of
this function should be reasonably
5580 *
template <
int dim>
5582 * ElastoPlasticProblem<dim>::solve_newton ()
5591 *
double residual_norm;
5592 *
double previous_residual_norm = -std::numeric_limits<double>::max();
5595 * previous_disp_norm = 0;
5597 *
const double correct_sigma = sigma_0;
5599 *
const unsigned int max_newton_iter = 100;
5601 *
for (
unsigned int newton_step = 1; newton_step <= max_newton_iter; ++newton_step)
5603 *
if (newton_step == 1
5605 * ((transfer_solution && timestep_no == 1)
5607 * !transfer_solution))
5608 * constitutive_law.set_sigma_0(1e+10);
5610 * constitutive_law.set_sigma_0(correct_sigma);
5612 * pcout <<
" " << std::endl;
5613 * pcout <<
" Newton iteration " << newton_step << std::endl;
5615 * pcout <<
" Assembling system... " << std::endl;
5616 * newton_matrix = 0;
5618 * newton_rhs_residual = 0;
5620 * tmp_solution = solution;
5621 * tmp_solution += incremental_displacement;
5622 * assemble_newton_system(tmp_solution,
5623 * incremental_displacement);
5625 * pcout <<
" Solving system... " << std::endl;
5626 * solve_newton_system();
5630 * It gets a bit more hairy after we have computed the
5631 * trial solution @f$\tilde{\mathbf u}@f$ of the current Newton step.
5632 * We handle a highly nonlinear problem so we have to damp
5633 * Newton
's method using a line search. To understand how we do this,
5634 * recall that in our formulation, we compute a trial solution
5635 * in each Newton step and not the update between old and new solution.
5636 * Since the solution set is a convex set, we will use a line
5637 * search that tries linear combinations of the
5638 * previous and the trial solution to guarantee that the
5639 * damped solution is in our solution set again.
5640 * At most we apply 5 damping steps.
5644 * There are exceptions to when we use a line search. First,
5645 * if this is the first Newton step on any mesh, then we don't have
5646 * any
point to compare the residual to, so we
always accept a full
5647 * step. Likewise,
if this is the
second Newton step on the
first mesh (or
5648 * the
second on any mesh
if we don
't transfer solutions from
5649 * mesh to mesh), then we have computed the first of these steps using
5650 * just an elastic model (see how we set the yield stress sigma to
5651 * an unreasonably large value above). In this case, the first Newton
5652 * solution was a purely elastic one, the second one a plastic one,
5653 * and any linear combination would not necessarily be expected to
5654 * lie in the feasible set -- so we just accept the solution we just
5659 * In either of these two cases, we bypass the line search and just
5660 * update residual and other vectors as necessary.
5663 * if ((newton_step==1)
5665 * (transfer_solution && newton_step == 2 && current_refinement_cycle == 0)
5667 * (!transfer_solution && newton_step == 2))
5669 * tmp_solution = solution;
5670 * tmp_solution += incremental_displacement;
5671 * compute_nonlinear_residual(tmp_solution);
5672 * old_solution = incremental_displacement;
5674 * residual = newton_rhs_residual;
5676 * residual.compress(VectorOperation::insert);
5678 * residual_norm = residual.l2_norm();
5680 * pcout << " Accepting Newton solution with residual: "
5681 * << residual_norm << std::endl;
5685 * for (unsigned int i = 0; i < 5; ++i)
5687 * distributed_solution = incremental_displacement;
5689 * const double alpha = std::pow(0.5, static_cast<double>(i));
5690 * tmp_vector = old_solution;
5691 * tmp_vector.sadd(1 - alpha, alpha, distributed_solution);
5693 * TimerOutput::Scope t(computing_timer, "Residual and lambda");
5695 * locally_relevant_tmp_vector = tmp_vector;
5696 * tmp_solution = solution;
5697 * tmp_solution += locally_relevant_tmp_vector;
5698 * compute_nonlinear_residual(tmp_solution);
5699 * residual = newton_rhs_residual;
5701 * residual.compress(VectorOperation::insert);
5703 * residual_norm = residual.l2_norm();
5705 * pcout << " Residual of the system: "
5706 * << residual_norm << std::endl
5707 * << " with a damping parameter alpha = " << alpha
5710 * if (residual_norm < previous_residual_norm)
5714 * incremental_displacement = tmp_vector;
5715 * old_solution = incremental_displacement;
5718 * disp_norm = incremental_displacement.l2_norm();
5723 * The final step is to check for convergence. If the residual is
5724 * less than a threshold of @f$10^{-10}@f$, then we terminate
5725 * the iteration on the current mesh:
5726 * if (residual_norm < 1e-10)
5729 * if (residual_norm < 1e-7)
5732 * pcout << " difference of two consecutive incremental displacement l2 norm : "
5733 * << std::abs(disp_norm - previous_disp_norm) << std::endl;
5734 * if ( std::abs(disp_norm - previous_disp_norm) < 1e-10 &&
5735 * (residual_norm < 1e-5 || std::abs(residual_norm - previous_residual_norm)<1e-9) )
5737 * pcout << " Convergence by difference of two consecutive solution! " << std::endl;
5742 * previous_residual_norm = residual_norm;
5743 * previous_disp_norm = disp_norm;
5750 * <a name="PlasticityContactProblemcompute_error"></a>
5751 * <h4>PlasticityContactProblem::compute_error</h4>
5757 * template <int dim>
5759 * ElastoPlasticProblem<dim>::compute_error ()
5761 * TrilinosWrappers::MPI::Vector tmp_solution(locally_owned_dofs, mpi_communicator);
5762 * tmp_solution = solution;
5763 * tmp_solution += incremental_displacement;
5765 * estimated_error_per_cell.reinit (triangulation.n_active_cells());
5766 * if (error_estimation_strategy == ErrorEstimationStrategy::kelly_error)
5768 * using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
5770 * KellyErrorEstimator<dim>::estimate(dof_handler,
5771 * QGauss<dim - 1>(fe.degree + 2),
5772 * std::map<types::boundary_id, const Function<dim> *>(),
5774 * estimated_error_per_cell);
5777 * else if (error_estimation_strategy == ErrorEstimationStrategy::residual_error)
5779 * compute_error_residual(tmp_solution);
5782 * else if (error_estimation_strategy == ErrorEstimationStrategy::weighted_residual_error)
5786 * make a non-parallel copy of tmp_solution
5789 * Vector<double> copy_solution(tmp_solution);
5793 * the dual function definition (it should be defined previously, e.g. input file)
5796 * if (base_mesh == "Timoshenko beam")
5798 * double length = .48,
5801 * const Point<dim> evaluation_point(length, -depth/2);
5803 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5805 * DualSolver<dim> dual_solver(triangulation, fe,
5807 * constitutive_law, dual_functional,
5808 * timestep_no, output_dir, base_mesh,
5809 * present_time, end_time);
5811 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5814 * else if (base_mesh == "Thick_tube_internal_pressure")
5816 * const unsigned int face_id = 0;
5817 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5818 * for (unsigned int i=0; i!=dim; ++i)
5820 * comp_stress[i].resize(dim);
5821 * for (unsigned int j=0; j!=dim; ++j)
5823 * comp_stress[i][j] = 1;
5827 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5829 * DualSolver<dim> dual_solver(triangulation, fe,
5831 * constitutive_law, dual_functional,
5832 * timestep_no, output_dir, base_mesh,
5833 * present_time, end_time);
5835 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5838 * else if (base_mesh == "Perforated_strip_tension")
5842 * .........................................
5843 * Mean stress_yy over the bottom boundary
5846 * const unsigned int face_id = 1;
5847 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5848 * for (unsigned int i=0; i!=dim; ++i)
5850 * comp_stress[i].resize(dim);
5851 * for (unsigned int j=0; j!=dim; ++j)
5853 * comp_stress[i][j] = 0;
5856 * comp_stress[1][1] = 1;
5858 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5862 * .........................................
5868 * DualSolver<dim> dual_solver(triangulation, fe,
5870 * constitutive_law, dual_functional,
5871 * timestep_no, output_dir, base_mesh,
5872 * present_time, end_time);
5874 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5877 * else if (base_mesh == "Cantiliver_beam_3d")
5881 * Quantity of interest:
5882 * -----------------------------------------------------------
5883 * displacement at Point A (x=0, y=height/2, z=length)
5887 * const double length = .7,
5890 * const Point<dim> evaluation_point(0, height/2, length);
5892 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5897 * -----------------------------------------------------------
5898 * Mean stress at the specified domain is of interest.
5899 * The interest domains are located on the bottom and top of the flanges
5900 * close to the clamped face, z = 0
5901 * top domain: height/2 - thickness_flange <= y <= height/2
5902 * 0 <= z <= 2 * thickness_flange
5903 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
5904 * 0 <= z <= 2 * thickness_flange
5910 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5911 * for (unsigned int i=0; i!=dim; ++i)
5913 * comp_stress[i].resize(dim);
5914 * for (unsigned int j=0; j!=dim; ++j)
5916 * comp_stress[i][j] = 1;
5919 * DualFunctional::MeanStressDomain<dim> dual_functional(base_mesh, comp_stress);
5923 * -----------------------------------------------------------
5929 * DualSolver<dim> dual_solver(triangulation, fe,
5931 * constitutive_law, dual_functional,
5932 * timestep_no, output_dir, base_mesh,
5933 * present_time, end_time);
5935 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5940 * AssertThrow(false, ExcNotImplemented());
5947 * AssertThrow(false, ExcNotImplemented());
5951 * relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm();
5953 * pcout << "Estimated relative error = " << relative_error << std::endl;
5957 * template <int dim>
5959 * ElastoPlasticProblem<dim>::compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution)
5961 * FEValues<dim> fe_values(fe, quadrature_formula,
5963 * update_gradients |
5965 * update_quadrature_points |
5966 * update_JxW_values);
5968 * const unsigned int n_q_points = quadrature_formula.size();
5969 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5970 * SymmetricTensor<4, dim> stress_strain_tensor_linearized;
5971 * SymmetricTensor<4, dim> stress_strain_tensor;
5972 * Tensor<5, dim> stress_strain_tensor_grad;
5973 * std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
5974 * for (unsigned int i=0; i!=n_q_points; ++i)
5976 * cell_hessians[i].resize (dim);
5978 * const EquationData::BodyForce<dim> body_force;
5980 * std::vector<Vector<double> > body_force_values (n_q_points, Vector<double>(dim));
5981 * const FEValuesExtractors::Vector displacement(0);
5984 * FEFaceValues<dim> fe_face_values_cell(fe, face_quadrature_formula,
5986 * update_quadrature_points|
5987 * update_gradients |
5988 * update_JxW_values |
5989 * update_normal_vectors),
5990 * fe_face_values_neighbor (fe, face_quadrature_formula,
5992 * update_gradients |
5993 * update_JxW_values |
5994 * update_normal_vectors);
5995 * FESubfaceValues<dim> fe_subface_values_cell (fe, face_quadrature_formula,
5996 * update_gradients);
5998 * const unsigned int n_face_q_points = face_quadrature_formula.size();
5999 * std::vector<Vector<double> > jump_residual (n_face_q_points, Vector<double>(dim));
6000 * std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
6001 * for (unsigned int i=0; i!=n_face_q_points; ++i)
6003 * cell_grads[i].resize (dim);
6005 * std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
6006 * for (unsigned int i=0; i!=n_face_q_points; ++i)
6008 * neighbor_grads[i].resize (dim);
6010 * SymmetricTensor<2, dim> q_cell_strain_tensor;
6011 * SymmetricTensor<2, dim> q_neighbor_strain_tensor;
6012 * SymmetricTensor<4, dim> cell_stress_strain_tensor;
6013 * SymmetricTensor<4, dim> neighbor_stress_strain_tensor;
6016 * typename std::map<typename DoFHandler<dim>::face_iterator, Vector<double> >
6018 * typename DoFHandler<dim>::active_cell_iterator
6019 * cell = dof_handler.begin_active(),
6020 * endc = dof_handler.end();
6021 * for (; cell!=endc; ++cell)
6022 * if (cell->is_locally_owned())
6024 * for (unsigned int face_no=0;
6025 * face_no<GeometryInfo<dim>::faces_per_cell;
6028 * face_integrals[cell->face(face_no)].reinit (dim);
6029 * face_integrals[cell->face(face_no)] = -1e20;
6033 * std::vector<Vector<float> > error_indicators_vector;
6034 * error_indicators_vector.resize( triangulation.n_active_cells(),
6035 * Vector<float>(dim) );
6039 * ----------------- estimate_some -------------------------
6042 * cell = dof_handler.begin_active();
6043 * unsigned int present_cell = 0;
6044 * for (; cell!=endc; ++cell, ++present_cell)
6045 * if (cell->is_locally_owned())
6049 * --------------- integrate_over_cell -------------------
6052 * fe_values.reinit(cell);
6053 * body_force.vector_value_list(fe_values.get_quadrature_points(),
6054 * body_force_values);
6055 * fe_values[displacement].get_function_symmetric_gradients(tmp_solution,
6057 * fe_values.get_function_hessians(tmp_solution, cell_hessians);
6059 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
6061 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
6062 * stress_strain_tensor_linearized,
6063 * stress_strain_tensor);
6064 * constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
6065 * cell_hessians[q_point],
6066 * stress_strain_tensor_grad);
6068 * for (unsigned int i=0; i!=dim; ++i)
6070 * error_indicators_vector[present_cell](i) +=
6071 * body_force_values[q_point](i)*fe_values.JxW(q_point);
6072 * for (unsigned int j=0; j!=dim; ++j)
6074 * for (unsigned int k=0; k!=dim; ++k)
6076 * for (unsigned int l=0; l!=dim; ++l)
6078 * error_indicators_vector[present_cell](i) +=
6079 * ( stress_strain_tensor[i][j][k][l]*
6080 * 0.5*(cell_hessians[q_point][k][l][j]
6082 * cell_hessians[q_point][l][k][j])
6083 * + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l]
6085 * fe_values.JxW(q_point);
6095 * -------------------------------------------------------
6096 * compute face_integrals
6099 * for (unsigned int face_no=0;
6100 * face_no<GeometryInfo<dim>::faces_per_cell;
6103 * if (cell->face(face_no)->at_boundary())
6105 * for (unsigned int id=0; id!=dim; ++id)
6107 * face_integrals[cell->face(face_no)](id) = 0;
6112 * if ((cell->neighbor(face_no)->has_children() == false) &&
6113 * (cell->neighbor(face_no)->level() == cell->level()) &&
6114 * (cell->neighbor(face_no)->index() < cell->index()))
6117 * if (cell->at_boundary(face_no) == false)
6118 * if (cell->neighbor(face_no)->level() < cell->level())
6122 * if (cell->face(face_no)->has_children() == false)
6126 * ------------- integrate_over_regular_face -----------
6129 * fe_face_values_cell.reinit(cell, face_no);
6130 * fe_face_values_cell.get_function_gradients (tmp_solution,
6133 * Assert (cell->neighbor(face_no).state() == IteratorState::valid,
6134 * ExcInternalError());
6135 * const unsigned int
6136 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6137 * const typename DoFHandler<dim>::active_cell_iterator
6138 * neighbor = cell->neighbor(face_no);
6140 * fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
6141 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6144 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6146 * q_cell_strain_tensor = 0.;
6147 * q_neighbor_strain_tensor = 0.;
6148 * for (unsigned int i=0; i!=dim; ++i)
6150 * for (unsigned int j=0; j!=dim; ++j)
6152 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6153 * cell_grads[q_point][j][i] );
6154 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6155 * neighbor_grads[q_point][j][i] );
6159 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6160 * cell_stress_strain_tensor);
6161 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6162 * neighbor_stress_strain_tensor);
6164 * jump_residual[q_point] = 0.;
6165 * for (unsigned int i=0; i!=dim; ++i)
6167 * for (unsigned int j=0; j!=dim; ++j)
6169 * for (unsigned int k=0; k!=dim; ++k)
6171 * for (unsigned int l=0; l!=dim; ++l)
6173 * jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
6174 * q_cell_strain_tensor[k][l]
6176 * neighbor_stress_strain_tensor[i][j][k][l]*
6177 * q_neighbor_strain_tensor[k][l] )*
6178 * fe_face_values_cell.normal_vector(q_point)[j];
6186 * Vector<double> face_integral_vector(dim);
6187 * face_integral_vector = 0;
6188 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6190 * for (unsigned int i=0; i!=dim; ++i)
6192 * face_integral_vector(i) += jump_residual[q_point](i) *
6193 * fe_face_values_cell.JxW(q_point);
6197 * Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
6198 * ExcInternalError());
6200 * for (unsigned int i=0; i!=dim; ++i)
6202 * Assert (face_integrals[cell->face(face_no)](i) == -1e20,
6203 * ExcInternalError());
6204 * face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
6210 * -----------------------------------------------------
6218 * ------------- integrate_over_irregular_face ---------
6221 * const typename DoFHandler<dim>::face_iterator
6222 * face = cell->face(face_no);
6223 * const typename DoFHandler<dim>::cell_iterator
6224 * neighbor = cell->neighbor(face_no);
6225 * Assert (neighbor.state() == IteratorState::valid,
6226 * ExcInternalError());
6227 * Assert (neighbor->has_children(),
6228 * ExcInternalError());
6230 * const unsigned int
6231 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6233 * for (unsigned int subface_no=0;
6234 * subface_no<face->n_children(); ++subface_no)
6236 * const typename DoFHandler<dim>::active_cell_iterator
6237 * neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
6238 * Assert (neighbor_child->face(neighbor_neighbor) ==
6239 * cell->face(face_no)->child(subface_no),
6240 * ExcInternalError());
6242 * fe_subface_values_cell.reinit (cell, face_no, subface_no);
6243 * fe_subface_values_cell.get_function_gradients (tmp_solution,
6245 * fe_face_values_neighbor.reinit (neighbor_child,
6246 * neighbor_neighbor);
6247 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6250 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6252 * q_cell_strain_tensor = 0.;
6253 * q_neighbor_strain_tensor = 0.;
6254 * for (unsigned int i=0; i!=dim; ++i)
6256 * for (unsigned int j=0; j!=dim; ++j)
6258 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6259 * cell_grads[q_point][j][i] );
6260 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6261 * neighbor_grads[q_point][j][i] );
6265 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6266 * cell_stress_strain_tensor);
6267 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6268 * neighbor_stress_strain_tensor);
6270 * jump_residual[q_point] = 0.;
6271 * for (unsigned int i=0; i!=dim; ++i)
6273 * for (unsigned int j=0; j!=dim; ++j)
6275 * for (unsigned int k=0; k!=dim; ++k)
6277 * for (unsigned int l=0; l!=dim; ++l)
6279 * jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
6280 * q_cell_strain_tensor[k][l]
6282 * neighbor_stress_strain_tensor[i][j][k][l]*
6283 * q_neighbor_strain_tensor[k][l] )*
6284 * fe_face_values_neighbor.normal_vector(q_point)[j];
6292 * Vector<double> face_integral_vector(dim);
6293 * face_integral_vector = 0;
6294 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6296 * for (unsigned int i=0; i!=dim; ++i)
6298 * face_integral_vector(i) += jump_residual[q_point](i) *
6299 * fe_face_values_neighbor.JxW(q_point);
6303 * for (unsigned int i=0; i!=dim; ++i)
6305 * face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
6310 * Vector<double> sum (dim);
6312 * for (unsigned int subface_no=0;
6313 * subface_no<face->n_children(); ++subface_no)
6315 * Assert (face_integrals.find(face->child(subface_no)) !=
6316 * face_integrals.end(),
6317 * ExcInternalError());
6318 * for (unsigned int i=0; i!=dim; ++i)
6320 * Assert (face_integrals[face->child(subface_no)](i) != -1e20,
6321 * ExcInternalError());
6322 * sum(i) += face_integrals[face->child(subface_no)](i);
6325 * for (unsigned int i=0; i!=dim; ++i)
6327 * face_integrals[face](i) = sum(i);
6333 * -----------------------------------------------------
6343 * ----------------------------------------------------------
6350 * cell = dof_handler.begin_active();
6351 * for (; cell!=endc; ++cell, ++present_cell)
6352 * if (cell->is_locally_owned())
6354 * for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
6357 * Assert(face_integrals.find(cell->face(face_no)) !=
6358 * face_integrals.end(),
6359 * ExcInternalError());
6361 * for (unsigned int id=0; id!=dim; ++id)
6363 * error_indicators_vector[present_cell](id)
6364 * -= 0.5*face_integrals[cell->face(face_no)](id);
6369 * estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
6379 * <a name="PlasticityContactProblemrefine_grid"></a>
6380 * <h4>PlasticityContactProblem::refine_grid</h4>
6384 * If you've made it
this far into the deal.II tutorial, the following
6385 * function refining the mesh should not pose any challenges to you
6386 * any more. It refines the mesh, either globally or
using the Kelly
6387 * error estimator, and
if so asked also transfers the solution from
6388 * the previous to the next mesh. In the latter
case, we also need
6389 * to compute the active set and other quantities again,
for which we
6390 * need the information computed by <code>compute_nonlinear_residual()</code>.
6393 * template <int dim>
6395 * ElastoPlasticProblem<dim>::refine_grid ()
6399 * ---------------------------------------------------------------
6400 * Make a field variable
for history varibales to be able to
6401 * transfer the data to the quadrature points of the
new mesh
6406 * history_dof_handler.distribute_dofs (history_fe);
6407 * std::vector< std::vector< Vector<double> > >
6408 * history_stress_field (dim, std::vector<
Vector<double> >(dim)),
6409 * local_history_stress_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6410 * local_history_stress_fe_values (dim, std::vector<
Vector<double> >(dim));
6413 * std::vector< std::vector< Vector<double> > >
6414 * history_strain_field (dim, std::vector<
Vector<double> >(dim)),
6415 * local_history_strain_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6416 * local_history_strain_fe_values (dim, std::vector<
Vector<double> >(dim));
6418 *
for (
unsigned int i=0; i<dim; ++i)
6419 *
for (
unsigned int j=0; j<dim; ++j)
6421 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6422 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6423 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6425 * history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6426 * local_history_strain_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6427 * local_history_strain_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6430 * quadrature_formula.size());
6433 * quadrature_formula, quadrature_formula,
6434 * qpoint_to_dof_matrix);
6436 * cell = dof_handler.begin_active(),
6437 * endc = dof_handler.end(),
6438 * dg_cell = history_dof_handler.begin_active();
6439 *
for (; cell!=endc; ++cell, ++dg_cell)
6440 *
if (cell->is_locally_owned())
6442 * PointHistory<dim> *local_quadrature_points_history
6443 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6444 *
Assert (local_quadrature_points_history >=
6445 * &quadrature_point_history.front(),
6446 * ExcInternalError());
6447 *
Assert (local_quadrature_points_history <
6448 * &quadrature_point_history.back(),
6449 * ExcInternalError());
6450 *
for (
unsigned int i=0; i<dim; ++i)
6451 *
for (
unsigned int j=0; j<dim; ++j)
6453 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6455 * local_history_stress_values_at_qpoints[i][j](q)
6456 * = local_quadrature_points_history[q].old_stress[i][j];
6458 * local_history_strain_values_at_qpoints[i][j](q)
6459 * = local_quadrature_points_history[q].old_strain[i][j];
6461 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
6462 * local_history_stress_values_at_qpoints[i][j]);
6463 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
6464 * history_stress_field[i][j]);
6466 * qpoint_to_dof_matrix.vmult (local_history_strain_fe_values[i][j],
6467 * local_history_strain_values_at_qpoints[i][j]);
6468 * dg_cell->set_dof_values (local_history_strain_fe_values[i][j],
6469 * history_strain_field[i][j]);
6476 * ---------------------------------------------------------------
6480 *
if (refinement_strategy == RefinementStrategy::refine_global)
6485 *
if (cell->is_locally_owned())
6486 * cell->set_refine_flag ();
6490 *
const double refine_fraction_cells = .3,
6491 * coarsen_fraction_cells = .03;
6494 *
const double refine_fraction_cells = .1,
6495 * coarsen_fraction_cells = .3;
6503 * estimated_error_per_cell,
6504 * refine_fraction_cells, coarsen_fraction_cells);
6511 * solution_transfer.prepare_for_coarsening_and_refinement(solution);
6516 *
if (transfer_solution)
6517 * incremental_displacement_transfer.prepare_for_coarsening_and_refinement(incremental_displacement);
6520 * history_stress_field_transfer1(history_dof_handler),
6521 * history_stress_field_transfer2(history_dof_handler);
6522 * history_stress_field_transfer0.prepare_for_coarsening_and_refinement(history_stress_field[0]);
6525 * history_stress_field_transfer1.prepare_for_coarsening_and_refinement(history_stress_field[1]);
6529 * history_stress_field_transfer2.prepare_for_coarsening_and_refinement(history_stress_field[2]);
6533 * history_strain_field_transfer1(history_dof_handler),
6534 * history_strain_field_transfer2(history_dof_handler);
6535 * history_strain_field_transfer0.prepare_for_coarsening_and_refinement(history_strain_field[0]);
6538 * history_strain_field_transfer1.prepare_for_coarsening_and_refinement(history_strain_field[1]);
6542 * history_strain_field_transfer2.prepare_for_coarsening_and_refinement(history_strain_field[2]);
6546 * pcout <<
" Number of active cells: "
6551 * setup_quadrature_point_history ();
6557 * distributed_solution = solution;
6560 * solution_transfer.interpolate(distributed_solution);
6561 * solution = distributed_solution;
6563 *
if (transfer_solution)
6568 * distributed_incremental_displacement = incremental_displacement;
6571 * incremental_displacement_transfer.interpolate(distributed_incremental_displacement);
6572 * incremental_displacement = distributed_incremental_displacement;
6575 * compute_nonlinear_residual(incremental_displacement);
6582 * ---------------------------------------------------
6585 * history_dof_handler.distribute_dofs (history_fe);
6591 * std::vector< std::vector< Vector<double> > >
6592 * distributed_history_stress_field (dim, std::vector<
Vector<double> >(dim));
6593 *
for (
unsigned int i=0; i<dim; ++i)
6594 *
for (
unsigned int j=0; j<dim; ++j)
6596 * distributed_history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6599 * history_stress_field_transfer0.interpolate(history_stress_field[0], distributed_history_stress_field[0]);
6602 * history_stress_field_transfer1.interpolate(history_stress_field[1], distributed_history_stress_field[1]);
6606 * history_stress_field_transfer2.interpolate(history_stress_field[2], distributed_history_stress_field[2]);
6609 * history_stress_field = distributed_history_stress_field;
6616 * std::vector< std::vector< Vector<double> > >
6617 * distributed_history_strain_field (dim, std::vector<
Vector<double> >(dim));
6618 *
for (
unsigned int i=0; i<dim; ++i)
6619 *
for (
unsigned int j=0; j<dim; ++j)
6621 * distributed_history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6624 * history_strain_field_transfer0.interpolate(history_strain_field[0], distributed_history_strain_field[0]);
6627 * history_strain_field_transfer1.interpolate(history_strain_field[1], distributed_history_strain_field[1]);
6631 * history_strain_field_transfer2.interpolate(history_strain_field[2], distributed_history_strain_field[2]);
6634 * history_strain_field = distributed_history_strain_field;
6638 * ---------------------------------------------------------------
6639 * Transfer the history data to the quadrature points of the
new mesh
6640 * In a
final step, we have to get the data back from the now
6641 * interpolated global field to the quadrature points on the
6642 *
new mesh. The following code will
do that:
6649 * history_fe.dofs_per_cell);
6652 * quadrature_formula,
6653 * dof_to_qpoint_matrix);
6654 * cell = dof_handler.begin_active();
6655 * endc = dof_handler.end();
6656 * dg_cell = history_dof_handler.begin_active();
6657 *
for (; cell != endc; ++cell, ++dg_cell)
6658 *
if (cell->is_locally_owned())
6660 * PointHistory<dim> *local_quadrature_points_history
6661 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6662 *
Assert (local_quadrature_points_history >=
6663 * &quadrature_point_history.front(),
6664 * ExcInternalError());
6665 *
Assert (local_quadrature_points_history <
6666 * &quadrature_point_history.back(),
6667 * ExcInternalError());
6668 *
for (
unsigned int i=0; i<dim; ++i)
6669 *
for (
unsigned int j=0; j<dim; ++j)
6671 * dg_cell->get_dof_values (history_stress_field[i][j],
6672 * local_history_stress_fe_values[i][j]);
6673 * dof_to_qpoint_matrix.vmult (local_history_stress_values_at_qpoints[i][j],
6674 * local_history_stress_fe_values[i][j]);
6676 * dg_cell->get_dof_values (history_strain_field[i][j],
6677 * local_history_strain_fe_values[i][j]);
6678 * dof_to_qpoint_matrix.vmult (local_history_strain_values_at_qpoints[i][j],
6679 * local_history_strain_fe_values[i][j]);
6680 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6682 * local_quadrature_points_history[q].old_stress[i][j]
6683 * = local_history_stress_values_at_qpoints[i][j](q);
6685 * local_quadrature_points_history[q].old_strain[i][j]
6686 * = local_history_strain_values_at_qpoints[i][j](q);
6697 * <a name=
"ElastoPlasticProblemsetup_quadrature_point_history"></a>
6698 * <h4>ElastoPlasticProblem::setup_quadrature_point_history</h4>
6703 * of the history variables, such as the existing stresses in the material,
6704 * that we store in each quadrature
point. As mentioned above, we use the
6705 * <code>user_pointer</code>
for this that is available in each cell.
6709 * To put
this into larger perspective, we note that
if we had previously
6710 * available stresses in our model (which we assume
do not exist
for the
6711 * purpose of
this program), then we would need to
interpolate the field of
6712 * preexisting stresses to the quadrature points. Likewise,
if we were to
6713 * simulate elasto-plastic materials with hardening/softening, then we would
6714 * have to store additional history variables like the present yield stress
6715 * of the accumulated plastic strains in each quadrature
6716 * points. Pre-existing hardening or weakening would then be implemented by
6717 * interpolating these variables in the present function as well.
6720 *
template <
int dim>
6721 *
void ElastoPlasticProblem<dim>::setup_quadrature_point_history ()
6725 * What we need to
do here is to
first count how many quadrature points
6726 * are within the responsibility of
this processor. This, of course,
6727 * equals the number of cells that belong to
this processor times the
6728 * number of quadrature points our quadrature formula has on each cell.
6732 * For good measure, we also
set all user pointers of all cells, whether
6733 * ours of not, to the null pointer. This way,
if we ever access the user
6734 * pointer of a cell which we should not have accessed, a segmentation
6735 * fault will let us know that
this should not have happened:
6738 *
unsigned int our_cells = 0;
6742 *
if (cell->is_locally_owned())
6749 * Next, allocate as many quadrature objects as we need. Since the
6750 * <code>resize</code> function does not actually shrink the amount of
6751 * allocated memory
if the requested
new size is smaller than the old
6752 * size, we resort to a trick to
first free all memory, and then
6753 * reallocate it: we declare an empty vector as a temporary variable and
6754 * then
swap the contents of the old vector and
this temporary
6755 * variable. This makes sure that the
6756 * <code>quadrature_point_history</code> is now really empty, and we can
6757 * let the temporary variable that now holds the previous contents of the
6758 * vector go out of scope and be destroyed. In the next step. we can then
6759 * re-allocate as many elements as we need, with the vector
6760 *
default-initializing the <code>PointHistory</code> objects, which
6761 * includes setting the stress variables to zero.
6765 * std::vector<PointHistory<dim> > tmp;
6766 * tmp.swap (quadrature_point_history);
6768 * quadrature_point_history.resize (our_cells *
6769 * quadrature_formula.size());
6773 * Finally
loop over all cells again and
set the user pointers from the
6774 * cells that belong to the present processor to
point to the
first
6775 * quadrature
point objects corresponding to
this cell in the vector of
6779 *
unsigned int history_index = 0;
6783 *
if (cell->is_locally_owned())
6785 * cell->set_user_pointer (&quadrature_point_history[history_index]);
6786 * history_index += quadrature_formula.size();
6791 * At the
end,
for good measure make sure that our count of elements was
6792 * correct and that we have both used up all objects we allocated
6793 * previously, and not
point to any objects beyond the
end of the
6794 * vector. Such defensive programming strategies are
always good checks to
6795 * avoid accidental errors and to guard against future changes to
this
6796 * function that forget to update all uses of a variable at the same
6797 * time. Recall that constructs
using the <code>
Assert</code> macro are
6798 * optimized away in optimized mode, so
do not affect the
run time of
6802 *
Assert (history_index == quadrature_point_history.size(),
6803 * ExcInternalError());
6809 * <a name=
"ElastoPlasticProblemupdate_quadrature_point_history"></a>
6810 * <h4>ElastoPlasticProblem::update_quadrature_point_history</h4>
6814 * At the
end of each time step, we should have computed an incremental
6815 * displacement update so that the material in its
new configuration
6816 * accommodates
for the difference between the external body and boundary
6817 * forces applied during
this time step minus the forces exerted through
6818 * preexisting
internal stresses. In order to have the preexisting
6819 * stresses available at the next time step, we therefore have to update the
6820 * preexisting stresses with the stresses due to the incremental
6821 * displacement computed during the present time step. Ideally, the
6822 * resulting
sum of
internal stresses would exactly counter all external
6823 * forces. Indeed, a simple experiment can make sure that
this is so:
if we
6824 * choose boundary conditions and body forces to be time independent, then
6825 * the forcing terms (the sum of external forces and
internal stresses)
6826 * should be exactly zero. If you make
this experiment, you will realize
6827 * from the output of the
norm of the right hand side in each time step that
6828 *
this is almost the
case: it is not exactly zero, since in the
first time
6829 * step the incremental displacement and stress updates were computed
6830 * relative to the undeformed mesh, which was then deformed. In the
second
6831 * time step, we again compute displacement and stress updates, but
this
6832 * time in the deformed mesh -- there, the resulting updates are very small
6833 * but not quite zero. This can be iterated, and in each such iteration the
6834 * residual, i.e. the
norm of the right hand side vector, is reduced;
if one
6835 * makes
this little experiment, one realizes that the
norm of
this residual
6836 * decays exponentially with the number of iterations, and after an
initial
6837 * very rapid decline is reduced by roughly a factor of about 3.5 in each
6838 * iteration (
for one testcase I looked at, other testcases, and other
6839 *
numbers of unknowns change the factor, but not the exponential decay).
6843 * In a sense,
this can then be considered as a quasi-timestepping scheme to
6844 * resolve the nonlinear problem of solving large-deformation elasticity on
6845 * a mesh that is moved along in a Lagrangian manner.
6849 * Another complication is that the existing (old) stresses are defined on
6850 * the old mesh, which we will move around after updating the stresses. If
6851 *
this mesh update involves rotations of the cell, then we need to also
6852 *
rotate the updated stress, since it was computed relative to the
6853 * coordinate system of the old cell.
6857 * Thus, what we need is the following: on each cell which the present
6858 * processor owns, we need to
extract the old stress from the data stored
6859 * with each quadrature
point, compute the stress update, add the two
6860 * together, and then
rotate the result together with the incremental
6861 * rotation computed from the incremental displacement at the present
6862 * quadrature
point. We will detail these steps below:
6865 *
template <
int dim>
6866 *
void ElastoPlasticProblem<dim>::
6867 * update_quadrature_point_history ()
6871 * First,
set up an <code>
FEValues</code>
object by which we will evaluate
6872 * the displacements and the
gradients thereof at the
6873 * quadrature points, together with a vector that will hold
this
6881 *
const unsigned int n_q_points = quadrature_formula.size();
6883 * std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
6889 * Then
loop over all cells and
do the job in the cells that belong to our
6897 * cell = dof_handler.begin_active(),
6898 * endc = dof_handler.end();
6902 *
for (; cell != endc; ++cell)
6903 *
if (cell->is_locally_owned())
6907 * Next, get a pointer to the quadrature
point history data local to
6908 * the present cell, and, as a defensive measure, make sure that
6909 *
this pointer is within the bounds of the global array:
6912 * PointHistory<dim> *local_quadrature_points_history
6913 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6914 *
Assert (local_quadrature_points_history >=
6915 * &quadrature_point_history.front(),
6916 * ExcInternalError());
6917 *
Assert (local_quadrature_points_history <
6918 * &quadrature_point_history.back(),
6919 * ExcInternalError());
6923 * Then initialize the <code>
FEValues</code>
object on the present
6924 * cell, and
extract the strains of the displacement at the
6928 * fe_values.reinit (cell);
6929 * fe_values[displacement].get_function_symmetric_gradients(incremental_displacement,
6930 * incremental_strain_tensor);
6934 * Then
loop over the quadrature points of
this cell:
6937 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6939 * local_quadrature_points_history[q].old_strain +=
6940 * incremental_strain_tensor[q];
6942 * constitutive_law.get_stress_strain_tensor(local_quadrature_points_history[q].old_strain,
6943 * stress_strain_tensor);
6947 * The result of these operations is then written back into
6948 * the original place:
6951 * local_quadrature_points_history[q].old_stress
6952 * = stress_strain_tensor * local_quadrature_points_history[q].old_strain;
6954 * local_quadrature_points_history[q].point
6955 * = fe_values.get_quadrature_points ()[q];
6964 * <a name=
"PlasticityContactProblemmove_mesh"></a>
6965 * <h4>PlasticityContactProblem::move_mesh</h4>
6969 * The remaining three
functions before we get to <code>
run()</code>
6970 * have to
do with generating output. The following one is an attempt
6971 * at showing the deformed body in its deformed configuration. To
this
6972 *
end,
this function takes a displacement vector field and moves every
6973 * vertex of the (local part) of the mesh by the previously computed
6974 * displacement. We will
call this function with the current
6975 * displacement field before we generate graphical output, and we will
6976 *
call it again after generating graphical output with the
negative
6977 * displacement field to undo the changes to the mesh so made.
6981 * The function itself is pretty straightforward. All we have to
do
6982 * is keep track which
vertices we have already touched, as we
6983 * encounter the same
vertices multiple times as we
loop over cells.
6986 *
template <
int dim>
6988 * ElastoPlasticProblem<dim>::
6991 * std::vector<bool> vertex_touched(
triangulation.n_vertices(),
false);
6994 * dof_handler.begin_active();
6995 * cell != dof_handler.end(); ++cell)
6996 *
if (cell->is_locally_owned())
6997 *
for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
6998 *
if (vertex_touched[cell->vertex_index(v)] ==
false)
7000 * vertex_touched[cell->vertex_index(v)] =
true;
7003 *
for (
unsigned int d = 0;
d < dim; ++
d)
7004 * vertex_displacement[d] = displacement(cell->vertex_dof_index(v, d));
7006 * cell->vertex(v) += vertex_displacement;
7015 * <a name=
"PlasticityContactProblemoutput_results"></a>
7016 * <h4>PlasticityContactProblem::output_results</h4>
7020 * Next is the function we use to actually generate graphical output. The
7021 * function is a bit tedious, but not actually particularly complicated.
7022 * It moves the mesh at the top (and moves it back at the end), then
7023 * computes the contact forces along the contact surface. We can
do
7024 * so (as shown in the accompanying paper) by taking the untreated
7025 * residual vector and identifying which degrees of freedom
7026 * correspond to those with contact by asking whether they have an
7027 * inhomogeneous constraints associated with them. As
always, we need
7028 * to be mindful that we can only write into completely distributed
7029 * vectors (i.e., vectors without ghost elements) but that when we
7030 * want to generate output, we need vectors that
do indeed have
7031 * ghost entries
for all locally relevant degrees of freedom.
7034 *
template <
int dim>
7036 * ElastoPlasticProblem<dim>::output_results (
const std::string &filename_base)
7040 * pcout <<
" Writing graphical output... " << std::flush;
7044 *
const double magnified_factor = 3;
7045 * magnified_solution *= magnified_factor;
7047 * move_mesh(magnified_solution);
7054 *
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
7056 * data_out.add_data_vector(solution,
7057 * std::vector<std::string> (dim,
"displacement"),
7061 * std::vector<std::string> solution_names;
7066 * solution_names.push_back (
"displacement");
7069 * solution_names.push_back (
"x_displacement");
7070 * solution_names.push_back (
"y_displacement");
7073 * solution_names.push_back (
"x_displacement");
7074 * solution_names.push_back (
"y_displacement");
7075 * solution_names.push_back (
"z_displacement");
7081 * data_out.add_data_vector (solution, solution_names);
7086 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
7088 * data_out.add_data_vector(subdomain,
"subdomain");
7091 * data_out.add_data_vector(fraction_of_plastic_q_points_per_cell,
7092 *
"fraction_of_plastic_q_points");
7095 * data_out.build_patches();
7099 * In the remainder of the function, we generate one VTU file on
7100 * every processor, indexed by the subdomain
id of
this processor.
7101 * On the
first processor, we then also create a <code>.pvtu</code>
7102 * file that indexes <i>all</i> of the VTU files so that the entire
7103 *
set of output files can be read at once. These <code>.pvtu</code>
7104 * are used by Paraview to describe an entire
parallel computation
's
7105 * output files. We then do the same again for the competitor of
7106 * Paraview, the Visit visualization program, by creating a matching
7107 * <code>.visit</code> file.
7110 * const std::string filename =
7111 * (output_dir + filename_base + "-"
7112 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7114 * std::ofstream output_vtu((filename + ".vtu").c_str());
7115 * data_out.write_vtu(output_vtu);
7116 * pcout << output_dir + filename_base << ".pvtu" << std::endl;
7119 * if (this_mpi_process == 0)
7121 * std::vector<std::string> filenames;
7122 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7123 * filenames.push_back(filename_base + "-" +
7124 * Utilities::int_to_string(i, 4) +
7127 * std::ofstream pvtu_master_output((output_dir + filename_base + ".pvtu").c_str());
7128 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7130 * std::ofstream visit_master_output((output_dir + filename_base + ".visit").c_str());
7131 * data_out.write_pvtu_record(visit_master_output, filenames);
7135 * produce eps files for mesh illustration
7138 * std::ofstream output_eps((filename + ".eps").c_str());
7140 * grid_out.write_eps(triangulation, output_eps);
7145 * Extrapolate the stresses from Gauss point to the nodes
7148 * SymmetricTensor<2, dim> stress_at_qpoint;
7150 * FE_DGQ<dim> history_fe (1);
7151 * DoFHandler<dim> history_dof_handler (triangulation);
7152 * history_dof_handler.distribute_dofs (history_fe);
7153 * std::vector< std::vector< Vector<double> > >
7154 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7155 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7156 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7157 * for (unsigned int i=0; i<dim; ++i)
7158 * for (unsigned int j=0; j<dim; ++j)
7160 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7161 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7162 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7165 * Vector<double> VM_stress_field (history_dof_handler.n_dofs()),
7166 * local_VM_stress_values_at_qpoints (quadrature_formula.size()),
7167 * local_VM_stress_fe_values (history_fe.dofs_per_cell);
7169 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7170 * quadrature_formula.size());
7171 * FETools::compute_projection_from_quadrature_points_matrix
7173 * quadrature_formula, quadrature_formula,
7174 * qpoint_to_dof_matrix);
7176 * typename DoFHandler<dim>::active_cell_iterator
7177 * cell = dof_handler.begin_active(),
7178 * endc = dof_handler.end(),
7179 * dg_cell = history_dof_handler.begin_active();
7181 * const FEValuesExtractors::Vector displacement(0);
7183 * for (; cell!=endc; ++cell, ++dg_cell)
7184 * if (cell->is_locally_owned())
7186 * PointHistory<dim> *local_quadrature_points_history
7187 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7188 * Assert (local_quadrature_points_history >=
7189 * &quadrature_point_history.front(),
7190 * ExcInternalError());
7191 * Assert (local_quadrature_points_history <
7192 * &quadrature_point_history.back(),
7193 * ExcInternalError());
7197 * Then loop over the quadrature points of this cell:
7200 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7202 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7204 * for (unsigned int i=0; i<dim; ++i)
7205 * for (unsigned int j=i; j<dim; ++j)
7207 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7210 * local_VM_stress_values_at_qpoints(q) = Evaluation::get_von_Mises_stress(stress_at_qpoint);
7215 * for (unsigned int i=0; i<dim; ++i)
7216 * for (unsigned int j=i; j<dim; ++j)
7218 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7219 * local_history_stress_values_at_qpoints[i][j]);
7220 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7221 * history_stress_field[i][j]);
7224 * qpoint_to_dof_matrix.vmult (local_VM_stress_fe_values,
7225 * local_VM_stress_values_at_qpoints);
7226 * dg_cell->set_dof_values (local_VM_stress_fe_values,
7234 * Save stresses on nodes by nodal averaging
7235 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7236 * in order to compute stresses on nodes (by applying nodal averaging)
7237 * Therefore, each vertex has one degree of freedom
7240 * FE_Q<dim> fe_1 (1);
7241 * DoFHandler<dim> dof_handler_1 (triangulation);
7242 * dof_handler_1.distribute_dofs (fe_1);
7244 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7245 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7247 * std::vector< std::vector< Vector<double> > >
7248 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7249 * for (unsigned int i=0; i<dim; ++i)
7250 * for (unsigned int j=0; j<dim; ++j)
7252 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7255 * Vector<double> VM_stress_on_vertices (dof_handler_1.n_dofs()),
7256 * counter_on_vertices (dof_handler_1.n_dofs());
7257 * VM_stress_on_vertices = 0;
7258 * counter_on_vertices = 0;
7260 * cell = dof_handler.begin_active();
7261 * dg_cell = history_dof_handler.begin_active();
7262 * typename DoFHandler<dim>::active_cell_iterator
7263 * cell_1 = dof_handler_1.begin_active();
7264 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7265 * if (cell->is_locally_owned())
7267 * dg_cell->get_dof_values (VM_stress_field,
7268 * local_VM_stress_fe_values);
7270 * for (unsigned int i=0; i<dim; ++i)
7271 * for (unsigned int j=0; j<dim; ++j)
7273 * dg_cell->get_dof_values (history_stress_field[i][j],
7274 * local_history_stress_fe_values[i][j]);
7277 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7279 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7284 * Point<dim> point1, point2;
7285 * point1 = cell_1->vertex(v);
7286 * point2 = dg_cell->vertex(v);
7287 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7294 * counter_on_vertices (dof_1_vertex) += 1;
7296 * VM_stress_on_vertices (dof_1_vertex) += local_VM_stress_fe_values (v);
7298 * for (unsigned int i=0; i<dim; ++i)
7299 * for (unsigned int j=0; j<dim; ++j)
7301 * history_stress_on_vertices[i][j](dof_1_vertex) +=
7302 * local_history_stress_fe_values[i][j](v);
7308 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7310 * VM_stress_on_vertices(id) /= counter_on_vertices(id);
7312 * for (unsigned int i=0; i<dim; ++i)
7313 * for (unsigned int j=0; j<dim; ++j)
7315 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7321 * Save figures of stresses
7324 * if (show_stresses)
7327 * DataOut<dim> data_out;
7328 * data_out.attach_dof_handler (history_dof_handler);
7331 * data_out.add_data_vector (history_stress_field[0][0], "stress_xx");
7332 * data_out.add_data_vector (history_stress_field[1][1], "stress_yy");
7333 * data_out.add_data_vector (history_stress_field[0][1], "stress_xy");
7334 * data_out.add_data_vector (VM_stress_field, "Von_Mises_stress");
7338 * data_out.add_data_vector (history_stress_field[0][2], "stress_xz");
7339 * data_out.add_data_vector (history_stress_field[1][2], "stress_yz");
7340 * data_out.add_data_vector (history_stress_field[2][2], "stress_zz");
7343 * data_out.build_patches ();
7345 * const std::string filename_base_stress = ("stress-" + filename_base);
7347 * const std::string filename =
7348 * (output_dir + filename_base_stress + "-"
7349 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7351 * std::ofstream output_vtu((filename + ".vtu").c_str());
7352 * data_out.write_vtu(output_vtu);
7353 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7355 * if (this_mpi_process == 0)
7357 * std::vector<std::string> filenames;
7358 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7359 * filenames.push_back(filename_base_stress + "-" +
7360 * Utilities::int_to_string(i, 4) +
7363 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7364 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7366 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7367 * data_out.write_pvtu_record(visit_master_output, filenames);
7374 * DataOut<dim> data_out;
7375 * data_out.attach_dof_handler (dof_handler_1);
7378 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_xx_averaged");
7379 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_yy_averaged");
7380 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_xy_averaged");
7381 * data_out.add_data_vector (VM_stress_on_vertices, "Von_Mises_stress_averaged");
7385 * data_out.add_data_vector (history_stress_on_vertices[0][2], "stress_xz_averaged");
7386 * data_out.add_data_vector (history_stress_on_vertices[1][2], "stress_yz_averaged");
7387 * data_out.add_data_vector (history_stress_on_vertices[2][2], "stress_zz_averaged");
7390 * data_out.build_patches ();
7392 * const std::string filename_base_stress = ("averaged-stress-" + filename_base);
7394 * const std::string filename =
7395 * (output_dir + filename_base_stress + "-"
7396 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7398 * std::ofstream output_vtu((filename + ".vtu").c_str());
7399 * data_out.write_vtu(output_vtu);
7400 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7402 * if (this_mpi_process == 0)
7404 * std::vector<std::string> filenames;
7405 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7406 * filenames.push_back(filename_base_stress + "-" +
7407 * Utilities::int_to_string(i, 4) +
7410 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7411 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7413 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7414 * data_out.write_pvtu_record(visit_master_output, filenames);
7421 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7429 * magnified_solution *= -1;
7430 * move_mesh(magnified_solution);
7437 * if (base_mesh == "Timoshenko beam")
7439 * const double length = .48,
7442 * Point<dim> intersted_point(length, -depth/2);
7443 * Point<dim> vertex_displacement;
7444 * bool vertex_found = false;
7446 * for (typename DoFHandler<dim>::active_cell_iterator cell =
7447 * dof_handler.begin_active();
7448 * cell != dof_handler.end(); ++cell)
7449 * if (cell->is_locally_owned() && !vertex_found)
7450 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7451 * if ( std::fabs(cell->vertex(v)[0] - intersted_point[0])<1e-6 &&
7452 * std::fabs(cell->vertex(v)[1] - intersted_point[1])<1e-6)
7454 * vertex_found = true;
7456 * for (unsigned int d = 0; d < dim; ++d)
7457 * vertex_displacement[d] = solution(cell->vertex_dof_index(v, d));
7462 * pcout << " Number of active cells: "
7463 * << triangulation.n_global_active_cells() << std::endl
7464 * << " Number of degrees of freedom: " << dof_handler.n_dofs()
7467 * AssertThrow(vertex_found, ExcInternalError());
7468 * std::cout << "Displacement at the point (" << intersted_point[0]
7469 * << ", " << intersted_point[1] << ") is "
7470 * << "(" << vertex_displacement[0]
7471 * << ", " << vertex_displacement[1] << ").\n";
7473 * Vector<double> vertex_exact_displacement(dim);
7474 * EquationData::IncrementalBoundaryValues<dim> incremental_boundary_values(present_time, end_time);
7475 * incremental_boundary_values.vector_value (intersted_point, vertex_exact_displacement);
7477 * std::cout << "Exact displacement at the point (" << intersted_point[0]
7478 * << ", " << intersted_point[1] << ") is "
7479 * << "(" << vertex_exact_displacement[0]
7480 * << ", " << vertex_exact_displacement[1] << ").\n\n";
7483 * else if (base_mesh == "Thick_tube_internal_pressure")
7485 * const double pressure (0.6*2.4e8),
7486 * inner_radius (.1);
7489 * const double pressure (1.94e8),
7490 * inner_radius (.1);
7498 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7499 * 3d and plane strain
7502 * const double mu (e_modulus / (2 * (1 + nu)));
7504 * const Point<dim> point_A(inner_radius, 0.);
7505 * Vector<double> disp_A(dim);
7509 * make a non-parallel copy of solution
7512 * Vector<double> copy_solution(solution);
7514 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
7516 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
7518 * table_results.add_value("time step", timestep_no);
7519 * table_results.add_value("Cells", triangulation.n_global_active_cells());
7520 * table_results.add_value("DoFs", dof_handler.n_dofs());
7521 * table_results.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7522 * table_results.add_value("4*mu*u_A/(sigma_0*a)", 4*mu*disp_A(0)/(sigma_0*inner_radius));
7526 * Compute stresses in the POLAR coordinates, 1- save it on Gauss points,
7527 * 2- extrapolate them to nodes and taking their avarages (nodal avaraging)
7530 * AssertThrow (dim == 2, ExcNotImplemented());
7534 * we define a rotation matrix to be able to transform the stress
7535 * from the Cartesian coordinate to the polar coordinate
7538 * Tensor<2, dim> rotation_matrix; // [cos sin; -sin cos] , sigma_r = rot * sigma * rot^T
7540 * FEValues<dim> fe_values (fe, quadrature_formula, update_quadrature_points |
7541 * update_values | update_gradients);
7543 * const unsigned int n_q_points = quadrature_formula.size();
7545 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
7546 * SymmetricTensor<4, dim> stress_strain_tensor;
7547 * Tensor<2, dim> stress_at_qpoint;
7549 * FE_DGQ<dim> history_fe (1);
7550 * DoFHandler<dim> history_dof_handler (triangulation);
7551 * history_dof_handler.distribute_dofs (history_fe);
7552 * std::vector< std::vector< Vector<double> > >
7553 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7554 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7555 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7556 * for (unsigned int i=0; i<dim; ++i)
7557 * for (unsigned int j=0; j<dim; ++j)
7559 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7560 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7561 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7564 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7565 * quadrature_formula.size());
7566 * FETools::compute_projection_from_quadrature_points_matrix
7568 * quadrature_formula, quadrature_formula,
7569 * qpoint_to_dof_matrix);
7571 * typename DoFHandler<dim>::active_cell_iterator
7572 * cell = dof_handler.begin_active(),
7573 * endc = dof_handler.end(),
7574 * dg_cell = history_dof_handler.begin_active();
7576 * const FEValuesExtractors::Vector displacement(0);
7578 * for (; cell!=endc; ++cell, ++dg_cell)
7579 * if (cell->is_locally_owned())
7581 * PointHistory<dim> *local_quadrature_points_history
7582 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7583 * Assert (local_quadrature_points_history >=
7584 * &quadrature_point_history.front(),
7585 * ExcInternalError());
7586 * Assert (local_quadrature_points_history <
7587 * &quadrature_point_history.back(),
7588 * ExcInternalError());
7592 * Then loop over the quadrature points of this cell:
7595 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7597 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7601 * transform the stress from the Cartesian coordinate to the polar coordinate
7604 * const Point<dim> point = local_quadrature_points_history[q].point;
7605 * const double theta = std::atan2(point(1),point(0));
7612 * rotation_matrix[0][0] = std::cos(theta);
7613 * rotation_matrix[0][1] = std::sin(theta);
7614 * rotation_matrix[1][0] = -std::sin(theta);
7615 * rotation_matrix[1][1] = std::cos(theta);
7619 * stress in polar coordinate
7622 * stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix);
7624 * for (unsigned int i=0; i<dim; ++i)
7625 * for (unsigned int j=i; j<dim; ++j)
7627 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7633 * for (unsigned int i=0; i<dim; ++i)
7634 * for (unsigned int j=i; j<dim; ++j)
7636 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7637 * local_history_stress_values_at_qpoints[i][j]);
7638 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7639 * history_stress_field[i][j]);
7645 * DataOut<dim> data_out;
7646 * data_out.attach_dof_handler (history_dof_handler);
7649 * data_out.add_data_vector (history_stress_field[0][0], "stress_rr");
7650 * data_out.add_data_vector (history_stress_field[1][1], "stress_tt");
7651 * data_out.add_data_vector (history_stress_field[0][1], "stress_rt");
7653 * data_out.build_patches ();
7655 * const std::string filename_base_stress = ("stress-polar-" + filename_base);
7657 * const std::string filename =
7658 * (output_dir + filename_base_stress + "-"
7659 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7661 * std::ofstream output_vtu((filename + ".vtu").c_str());
7662 * data_out.write_vtu(output_vtu);
7663 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7665 * if (this_mpi_process == 0)
7667 * std::vector<std::string> filenames;
7668 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7669 * filenames.push_back(filename_base_stress + "-" +
7670 * Utilities::int_to_string(i, 4) +
7673 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7674 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7676 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7677 * data_out.write_pvtu_record(visit_master_output, filenames);
7685 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7686 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7687 * in order to compute stresses on nodes (by applying nodal averaging)
7688 * Therefore, each vertex has one degree of freedom
7691 * FE_Q<dim> fe_1 (1);
7692 * DoFHandler<dim> dof_handler_1 (triangulation);
7693 * dof_handler_1.distribute_dofs (fe_1);
7695 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7696 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7698 * std::vector< std::vector< Vector<double> > >
7699 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7700 * for (unsigned int i=0; i<dim; ++i)
7701 * for (unsigned int j=0; j<dim; ++j)
7703 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7706 * Vector<double> counter_on_vertices (dof_handler_1.n_dofs());
7707 * counter_on_vertices = 0;
7709 * cell = dof_handler.begin_active();
7710 * dg_cell = history_dof_handler.begin_active();
7711 * typename DoFHandler<dim>::active_cell_iterator
7712 * cell_1 = dof_handler_1.begin_active();
7713 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7714 * if (cell->is_locally_owned())
7717 * for (unsigned int i=0; i<dim; ++i)
7718 * for (unsigned int j=0; j<dim; ++j)
7720 * dg_cell->get_dof_values (history_stress_field[i][j],
7721 * local_history_stress_fe_values[i][j]);
7724 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7726 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7731 * Point<dim> point1, point2;
7732 * point1 = cell_1->vertex(v);
7733 * point2 = dg_cell->vertex(v);
7734 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7741 * counter_on_vertices (dof_1_vertex) += 1;
7743 * for (unsigned int i=0; i<dim; ++i)
7744 * for (unsigned int j=0; j<dim; ++j)
7746 * history_stress_on_vertices[i][j](dof_1_vertex) +=
7747 * local_history_stress_fe_values[i][j](v);
7753 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7755 * for (unsigned int i=0; i<dim; ++i)
7756 * for (unsigned int j=0; j<dim; ++j)
7758 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7764 * DataOut<dim> data_out;
7765 * data_out.attach_dof_handler (dof_handler_1);
7768 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_rr_averaged");
7769 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_tt_averaged");
7770 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_rt_averaged");
7772 * data_out.build_patches ();
7774 * const std::string filename_base_stress = ("averaged-stress-polar-" + filename_base);
7776 * const std::string filename =
7777 * (output_dir + filename_base_stress + "-"
7778 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7780 * std::ofstream output_vtu((filename + ".vtu").c_str());
7781 * data_out.write_vtu(output_vtu);
7782 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7784 * if (this_mpi_process == 0)
7786 * std::vector<std::string> filenames;
7787 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7788 * filenames.push_back(filename_base_stress + "-" +
7789 * Utilities::int_to_string(i, 4) +
7792 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7793 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7795 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7796 * data_out.write_pvtu_record(visit_master_output, filenames);
7803 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7809 * if ( std::abs( (present_time/end_time)*(pressure/sigma_0) - 0.6 ) <
7810 * .501*(present_timestep/end_time)*(pressure/sigma_0) )
7815 * table_results_2: presenting the stress_rr and stress_tt on the nodes of bottom edge
7818 * const unsigned int face_id = 3;
7820 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
7822 * bool evaluation_face_found = false;
7824 * typename DoFHandler<dim>::active_cell_iterator
7825 * cell = dof_handler.begin_active(),
7826 * endc = dof_handler.end(),
7827 * cell_1 = dof_handler_1.begin_active();
7828 * for (; cell!=endc; ++cell, ++cell_1)
7829 * if (cell->is_locally_owned())
7831 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7833 * if (cell->face(face)->at_boundary()
7835 * cell->face(face)->boundary_id() == face_id)
7837 * if (!evaluation_face_found)
7839 * evaluation_face_found = true;
7843 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7845 * types::global_dof_index dof_1_vertex =
7846 * cell_1->face(face)->vertex_dof_index(v, 0);
7847 * if (!vertices_found[dof_1_vertex])
7850 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
7852 * table_results_2.add_value("x coordinate", vertex_coordinate[0]);
7853 * table_results_2.add_value("stress_rr", history_stress_on_vertices[0][0](dof_1_vertex));
7854 * table_results_2.add_value("stress_tt", history_stress_on_vertices[1][1](dof_1_vertex));
7855 * table_results_2.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7857 * vertices_found[dof_1_vertex] = true;
7866 * AssertThrow(evaluation_face_found, ExcInternalError());
7870 * table_results_3: presenting the mean stress_rr of the nodes on the inner radius
7873 * const unsigned int face_id_2 = 0;
7875 * Tensor<2, dim> stress_node,
7876 * mean_stress_polar;
7877 * mean_stress_polar = 0;
7879 * std::vector<bool> vertices_found_2 (dof_handler_1.n_dofs(), false);
7880 * unsigned int no_vertices_found = 0;
7882 * evaluation_face_found = false;
7884 * cell = dof_handler.begin_active(),
7885 * endc = dof_handler.end(),
7886 * cell_1 = dof_handler_1.begin_active();
7887 * for (; cell!=endc; ++cell, ++cell_1)
7888 * if (cell->is_locally_owned())
7890 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7892 * if (cell->face(face)->at_boundary()
7894 * cell->face(face)->boundary_id() == face_id_2)
7896 * if (!evaluation_face_found)
7898 * evaluation_face_found = true;
7902 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7904 * types::global_dof_index dof_1_vertex =
7905 * cell_1->face(face)->vertex_dof_index(v, 0);
7906 * if (!vertices_found_2[dof_1_vertex])
7908 * for (unsigned int ir=0; ir<dim; ++ir)
7909 * for (unsigned int ic=0; ic<dim; ++ic)
7910 * stress_node[ir][ic] = history_stress_on_vertices[ir][ic](dof_1_vertex);
7912 * mean_stress_polar += stress_node;
7914 * vertices_found_2[dof_1_vertex] = true;
7915 * ++no_vertices_found;
7924 * AssertThrow(evaluation_face_found, ExcInternalError());
7926 * mean_stress_polar /= no_vertices_found;
7928 * table_results_3.add_value("time step", timestep_no);
7929 * table_results_3.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7930 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
7931 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
7932 * table_results_3.add_value("radius", inner_radius);
7933 * table_results_3.add_value("mean stress_rr", mean_stress_polar[0][0]);
7934 * table_results_3.add_value("mean stress_tt", mean_stress_polar[1][1]);
7941 * else if (base_mesh == "Perforated_strip_tension")
7943 * const double imposed_displacement (0.00055),
7944 * inner_radius (0.05);
7949 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7950 * 3d and plane strain
7954 * table_results: Demonstrates the result of displacement at the top left corner versus imposed tension
7959 * const Point<dim> point_C(0., height);
7960 * Vector<double> disp_C(dim);
7964 * make a non-parallel copy of solution
7967 * Vector<double> copy_solution(solution);
7969 * typename Evaluation::PointValuesEvaluation<dim>::
7970 * PointValuesEvaluation point_values_evaluation(point_C);
7972 * point_values_evaluation.compute (dof_handler, copy_solution, disp_C);
7974 * table_results.add_value("time step", timestep_no);
7975 * table_results.add_value("Cells", triangulation.n_global_active_cells());
7976 * table_results.add_value("DoFs", dof_handler.n_dofs());
7977 * table_results.add_value("4*mu*u_C/(sigma_0*r)", 4*mu*disp_C(1)/(sigma_0*inner_radius));
7983 * compute average sigma_yy on the bottom edge
7986 * double stress_yy_av;
7989 * const unsigned int face_id = 1;
7991 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
7992 * unsigned int no_vertices_in_face = 0;
7994 * bool evaluation_face_found = false;
7996 * typename DoFHandler<dim>::active_cell_iterator
7997 * cell = dof_handler.begin_active(),
7998 * endc = dof_handler.end(),
7999 * cell_1 = dof_handler_1.begin_active();
8000 * for (; cell!=endc; ++cell, ++cell_1)
8001 * if (cell->is_locally_owned())
8003 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
8005 * if (cell->face(face)->at_boundary()
8007 * cell->face(face)->boundary_id() == face_id)
8009 * if (!evaluation_face_found)
8011 * evaluation_face_found = true;
8015 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
8017 * types::global_dof_index dof_1_vertex =
8018 * cell_1->face(face)->vertex_dof_index(v, 0);
8019 * if (!vertices_found[dof_1_vertex])
8021 * stress_yy_av += history_stress_on_vertices[1][1](dof_1_vertex);
8022 * ++no_vertices_in_face;
8024 * vertices_found[dof_1_vertex] = true;
8033 * AssertThrow(evaluation_face_found, ExcInternalError());
8035 * stress_yy_av /= no_vertices_in_face;
8041 * table_results_2: Demonstrate the stress_yy on the nodes of bottom edge
8045 * if ( std::abs( (stress_yy_av/sigma_0) - .91 ) < .2 )
8048 * if ( (timestep_no) % 19 == 0 )
8055 * const unsigned int face_id = 1;
8057 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
8059 * bool evaluation_face_found = false;
8061 * typename DoFHandler<dim>::active_cell_iterator
8062 * cell = dof_handler.begin_active(),
8063 * endc = dof_handler.end(),
8064 * cell_1 = dof_handler_1.begin_active();
8065 * for (; cell!=endc; ++cell, ++cell_1)
8066 * if (cell->is_locally_owned())
8068 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
8070 * if (cell->face(face)->at_boundary()
8072 * cell->face(face)->boundary_id() == face_id)
8074 * if (!evaluation_face_found)
8076 * evaluation_face_found = true;
8080 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
8082 * types::global_dof_index dof_1_vertex =
8083 * cell_1->face(face)->vertex_dof_index(v, 0);
8085 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
8087 * if (!vertices_found[dof_1_vertex] && std::abs(vertex_coordinate[2])<1.e-8)
8089 * table_results_2.add_value("x", vertex_coordinate[0]);
8090 * table_results_2.add_value("x/r", vertex_coordinate[0]/inner_radius);
8091 * table_results_2.add_value("stress_xx/sigma_0", history_stress_on_vertices[0][0](dof_1_vertex)/sigma_0);
8092 * table_results_2.add_value("stress_yy/sigma_0", history_stress_on_vertices[1][1](dof_1_vertex)/sigma_0);
8093 * table_results_2.add_value("stress_yy_av/sigma_0", stress_yy_av/sigma_0);
8094 * table_results_2.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8096 * vertices_found[dof_1_vertex] = true;
8105 * AssertThrow(evaluation_face_found, ExcInternalError());
8111 * table_results_3: Demonstrate the Stress_mean (average tensile stress)
8112 * on the bottom edge versus epsilon_yy on the bottom left corner
8116 * double strain_yy_A = 0.;
8120 * compute strain_yy_A
8121 * Since the point A is the node on the bottom left corner,
8122 * we need to work just with one element
8126 * const Point<dim> point_A(inner_radius, 0, 0);
8128 * Vector<double> local_strain_yy_values_at_qpoints (quadrature_formula.size()),
8129 * local_strain_yy_fe_values (history_fe.dofs_per_cell);
8131 * SymmetricTensor<2, dim> strain_at_qpoint;
8133 * typename DoFHandler<dim>::active_cell_iterator
8134 * cell = dof_handler.begin_active(),
8135 * endc = dof_handler.end(),
8136 * dg_cell = history_dof_handler.begin_active();
8138 * bool cell_found = false;
8140 * for (; cell!=endc; ++cell, ++dg_cell)
8141 * if (cell->is_locally_owned() && !cell_found)
8143 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
8144 * if ( std::fabs(cell->vertex(v)[0] - point_A[0])<1e-6 &&
8145 * std::fabs(cell->vertex(v)[1] - point_A[1])<1e-6 &&
8146 * std::fabs(cell->vertex(v)[2] - point_A[2])<1e-6)
8148 * PointHistory<dim> *local_quadrature_points_history
8149 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8150 * Assert (local_quadrature_points_history >=
8151 * &quadrature_point_history.front(),
8152 * ExcInternalError());
8153 * Assert (local_quadrature_points_history <
8154 * &quadrature_point_history.back(),
8155 * ExcInternalError());
8159 * Then loop over the quadrature points of this cell:
8162 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8164 * strain_at_qpoint = local_quadrature_points_history[q].old_strain;
8166 * local_strain_yy_values_at_qpoints(q) = strain_at_qpoint[1][1];
8169 * qpoint_to_dof_matrix.vmult (local_strain_yy_fe_values,
8170 * local_strain_yy_values_at_qpoints);
8172 * strain_yy_A = local_strain_yy_fe_values (v);
8174 * cell_found = true;
8182 * table_results_3.add_value("time step", timestep_no);
8183 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
8184 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
8185 * table_results_3.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8186 * table_results_3.add_value("mean_tensile_stress/sigma_0", stress_yy_av/sigma_0);
8187 * table_results_3.add_value("E*strain_yy-A/sigma_0", e_modulus*strain_yy_A/sigma_0);
8192 * if (std::abs(present_time-end_time) < 1.e-7)
8194 * table_results_2.set_precision("Imposed u_y", 6);
8195 * table_results_3.set_precision("Imposed u_y", 6);
8199 * else if (base_mesh == "Cantiliver_beam_3d")
8201 * const double pressure (6e6),
8207 * table_results: Demonstrates the result of displacement at the top front point, Point A
8213 * Quantity of interest:
8214 * displacement at Point A (x=0, y=height/2, z=length)
8220 * const Point<dim> point_A(0, height/2, length);
8221 * Vector<double> disp_A(dim);
8225 * make a non-parallel copy of solution
8228 * Vector<double> copy_solution(solution);
8230 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
8232 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
8234 * table_results.add_value("time step", timestep_no);
8235 * table_results.add_value("Cells", triangulation.n_global_active_cells());
8236 * table_results.add_value("DoFs", dof_handler.n_dofs());
8237 * table_results.add_value("pressure", pressure*present_time/end_time);
8238 * table_results.add_value("u_A", disp_A(1));
8244 * demonstrate the location and maximum von-Mises stress in the
8245 * specified domain close to the clamped face, z = 0
8246 * top domain: height/2 - thickness_flange <= y <= height/2
8247 * 0 <= z <= 2 * thickness_flange
8248 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
8249 * 0 <= z <= 2 * thickness_flange
8255 * double VM_stress_max (0);
8256 * Point<dim> point_max;
8258 * SymmetricTensor<2, dim> stress_at_qpoint;
8260 * typename DoFHandler<dim>::active_cell_iterator
8261 * cell = dof_handler.begin_active(),
8262 * endc = dof_handler.end();
8264 * const FEValuesExtractors::Vector displacement(0);
8266 * for (; cell!=endc; ++cell)
8267 * if (cell->is_locally_owned())
8269 * PointHistory<dim> *local_quadrature_points_history
8270 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8271 * Assert (local_quadrature_points_history >=
8272 * &quadrature_point_history.front(),
8273 * ExcInternalError());
8274 * Assert (local_quadrature_points_history <
8275 * &quadrature_point_history.back(),
8276 * ExcInternalError());
8280 * Then loop over the quadrature points of this cell:
8283 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8285 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
8287 * const double VM_stress = Evaluation::get_von_Mises_stress(stress_at_qpoint);
8288 * if (VM_stress > VM_stress_max)
8290 * VM_stress_max = VM_stress;
8291 * point_max = local_quadrature_points_history[q].point;
8297 * table_results.add_value("maximum von_Mises stress", VM_stress_max);
8298 * table_results.add_value("x", point_max[0]);
8299 * table_results.add_value("y", point_max[1]);
8300 * table_results.add_value("z", point_max[2]);
8313 * <a name="PlasticityContactProblemrun"></a>
8314 * <h4>PlasticityContactProblem::run</h4>
8318 * As in all other tutorial programs, the <code>run()</code> function contains
8319 * the overall logic. There is not very much to it here: in essence, it
8320 * performs the loops over all mesh refinement cycles, and within each, hands
8321 * things over to the Newton solver in <code>solve_newton()</code> on the
8322 * current mesh and calls the function that creates graphical output for
8323 * the so-computed solution. It then outputs some statistics concerning both
8324 * run times and memory consumption that has been collected over the course of
8325 * computations on this mesh.
8328 * template <int dim>
8330 * ElastoPlasticProblem<dim>::run ()
8332 * computing_timer.reset();
8335 * present_timestep = 1;
8343 * ----------------------------------------------------------------
8344 * base_mesh == "Thick_tube_internal_pressure"
8348 * const Point<dim> center(0, 0);
8349 * const double inner_radius = .1,
8350 * outer_radius = .2;
8352 * const SphericalManifold<dim> inner_boundary_description(center, inner_radius);
8353 * triangulation.set_manifold (0, inner_boundary_description);
8355 * const SphericalManifold<dim> outer_boundary_description(center, outer_radius);
8356 * triangulation.set_manifold (1, outer_boundary_description);
8360 * ----------------------------------------------------------------
8361 * base_mesh == "Perforated_strip_tension"
8365 * const double inner_radius = 0.05;
8367 * const CylinderBoundary<dim> inner_boundary_description(inner_radius, 2);
8368 * triangulation.set_manifold (10, inner_boundary_description);
8372 * ----------------------------------------------------------------
8378 * setup_quadrature_point_history ();
8380 * while (present_time < end_time)
8382 * present_time += present_timestep;
8385 * if (present_time > end_time)
8387 * present_timestep -= (present_time - end_time);
8388 * present_time = end_time;
8390 * pcout << std::endl;
8391 * pcout << "Time step " << timestep_no << " at time " << present_time
8394 * relative_error = max_relative_error * 10;
8395 * current_refinement_cycle = 0;
8402 * ------------------------ Refinement based on the relative error -------------------------------
8408 * while (relative_error >= max_relative_error)
8413 * if ( (timestep_no > 1) && (current_refinement_cycle>0) && (relative_error >= max_relative_error) )
8415 * pcout << "The relative error, " << relative_error
8416 * << " , is still more than maximum relative error, "
8417 * << max_relative_error << ", but we move to the next increment.\n";
8418 * relative_error = .1 * max_relative_error;
8421 * if (relative_error >= max_relative_error)
8423 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8424 * ++current_refinement_cycle;
8432 * ------------------------ Refinement based on the number of refinement --------------------------
8436 * bool continue_loop = true;
8437 * while (continue_loop)
8442 * if ( (timestep_no == 1) && (current_refinement_cycle < 1) )
8444 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8445 * ++current_refinement_cycle;
8449 * continue_loop = false;
8457 * -------------------------------------------------------------------------------------------------
8463 * solution += incremental_displacement;
8465 * update_quadrature_point_history ();
8467 * output_results((std::string("solution-") +
8468 * Utilities::int_to_string(timestep_no, 4)).c_str());
8470 * computing_timer.print_summary();
8471 * computing_timer.reset();
8473 * Utilities::System::MemoryStats stats;
8474 * Utilities::System::get_memory_stats(stats);
8475 * pcout << "Peak virtual memory used, resident in kB: " << stats.VmSize << " "
8476 * << stats.VmRSS << std::endl;
8479 * if (std::abs(present_time-end_time) < 1.e-7)
8481 * const std::string filename = (output_dir + "Results");
8483 * std::ofstream output_txt((filename + ".txt").c_str());
8485 * pcout << std::endl;
8486 * table_results.write_text(output_txt);
8487 * pcout << std::endl;
8488 * table_results_2.write_text(output_txt);
8489 * pcout << std::endl;
8490 * table_results_3.write_text(output_txt);
8491 * pcout << std::endl;
8496 * if (base_mesh == "Thick_tube_internal_pressure")
8498 * triangulation.reset_manifold (0);
8499 * triangulation.reset_manifold (1);
8501 * else if (base_mesh == "Perforated_strip_tension")
8503 * triangulation.reset_manifold (10);
8512 * <a name="Thecodemaincodefunction"></a>
8513 * <h3>The <code>main</code> function</h3>
8517 * There really isn't much to the <code>main()</code> function. It looks
8521 *
int main (
int argc,
char *argv[])
8523 *
using namespace dealii;
8524 *
using namespace ElastoPlastic;
8530 *
const int dim = 3;
8531 * ElastoPlasticProblem<dim>::declare_parameters(prm);
8534 * std::cerr <<
"*** Call this program as <./elastoplastic input.prm>" << std::endl;
8538 * prm.parse_input(argv[1]);
8541 * ElastoPlasticProblem<dim> problem(prm);
8545 *
catch (std::exception &exc)
8547 * std::cerr << std::endl << std::endl
8548 * <<
"----------------------------------------------------"
8550 * std::cerr <<
"Exception on processing: " << std::endl
8551 * << exc.what() << std::endl
8552 * <<
"Aborting!" << std::endl
8553 * <<
"----------------------------------------------------"
8560 * std::cerr << std::endl << std::endl
8561 * <<
"----------------------------------------------------"
8563 * std::cerr <<
"Unknown exception!" << std::endl
8564 * <<
"Aborting!" << std::endl
8565 * <<
"----------------------------------------------------"
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
__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(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
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=AffineConstraints< number >(), 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 apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
@ component_is_part_of_vector
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
SignedDistance::Plane< dim > Plane
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)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
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)
::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< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
const ::Triangulation< dim, spacedim > & tria
1.9.7