158 * <a name=
"elastoplastic.cc-Includefiles"></a>
159 * <h3>Include files</h3>
160 * The set of include files is not much of a surprise any more at
this time:
163 * #include <deal.II/base/conditional_ostream.h>
164 * #include <deal.II/base/parameter_handler.h>
165 * #include <deal.II/base/utilities.h>
166 * #include <deal.II/base/index_set.h>
167 * #include <deal.II/base/quadrature_lib.h>
168 * #include <deal.II/base/function.h>
169 * #include <deal.II/base/logstream.h>
170 * #include <deal.II/base/timer.h>
171 * #include <deal.II/base/table_handler.h>
173 * #include <deal.II/lac/vector.h>
174 * #include <deal.II/lac/full_matrix.h>
175 * #include <deal.II/lac/sparsity_tools.h>
176 * #include <deal.II/lac/sparse_matrix.h>
177 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
178 * #include <deal.II/lac/block_sparsity_pattern.h>
179 * #include <deal.II/lac/solver_bicgstab.h>
180 * #include <deal.II/lac/precondition.h>
181 * #include <deal.II/lac/affine_constraints.h>
182 * #include <deal.II/lac/trilinos_sparse_matrix.h>
183 * #include <deal.II/lac/trilinos_block_sparse_matrix.h>
184 * #include <deal.II/lac/trilinos_vector.h>
185 * #include <deal.II/lac/trilinos_precondition.h>
186 * #include <deal.II/lac/trilinos_solver.h>
187 * #include <deal.II/lac/sparse_direct.h>
189 * #include <deal.II/grid/tria.h>
190 * #include <deal.II/grid/grid_generator.h>
191 * #include <deal.II/grid/grid_refinement.h>
192 * #include <deal.II/grid/grid_tools.h>
193 * #include <deal.II/grid/tria_accessor.h>
194 * #include <deal.II/grid/tria_iterator.h>
195 * #include <deal.II/grid/grid_out.h>
196 * #include <deal.II/grid/manifold_lib.h>
198 * #include <deal.II/distributed/tria.h>
199 * #include <deal.II/distributed/grid_refinement.h>
200 * #include <deal.II/distributed/solution_transfer.h>
202 * #include <deal.II/dofs/dof_handler.h>
203 * #include <deal.II/dofs/dof_accessor.h>
204 * #include <deal.II/dofs/dof_renumbering.h>
205 * #include <deal.II/dofs/dof_tools.h>
207 * #include <deal.II/fe/fe_q.h>
208 * #include <deal.II/fe/fe_system.h>
209 * #include <deal.II/fe/fe_values.h>
210 * #include <deal.II/fe/fe_dgq.h>
211 * #include <deal.II/fe/fe_tools.h>
213 * #include <deal.II/numerics/vector_tools.h>
214 * #include <deal.II/numerics/matrix_tools.h>
215 * #include <deal.II/numerics/data_out.h>
216 * #include <deal.II/numerics/error_estimator.h>
217 * #include <deal.II/numerics/fe_field_function.h>
218 * #include <deal.II/numerics/solution_transfer.h>
222 * And here the only two
new things among the header files: an include file in
223 * which
symmetric tensors of rank 2 and 4 are implemented, as introduced in
227 * #include <deal.II/base/symmetric_tensor.h>
231 * And a header that implements filters
for iterators looping over all
232 * cells. We will use
this when selecting only those cells
for output that are
233 * owned by the present process in a %
parallel program:
236 * #include <deal.II/grid/filtered_iterator.h>
239 * #include <iostream>
243 * This
final include file provides the <code>mkdir</code> function
244 * that we will use to create a directory
for output files,
if necessary:
247 * #include <sys/stat.h>
249 *
namespace ElastoPlastic
255 *
const unsigned int n_slices,
256 *
const double height,
261 *
Assert (input.n_levels() == 1,
262 * ExcMessage (
"The input triangulations must be coarse meshes."));
265 *
Assert(result.n_cells()==0, ExcMessage(
"resultin Triangulation need to be empty upon calling extrude_triangulation."));
266 *
Assert(height>0, ExcMessage(
"The height in extrude_triangulation needs to be positive."));
267 *
Assert(n_slices>=2, ExcMessage(
"The number of slices in extrude_triangulation needs to be at least 2."));
269 * std::vector<Point<3> > points(n_slices*input.n_used_vertices());
270 * std::vector<CellData<3> > cells;
271 * cells.reserve((n_slices-1)*input.n_active_cells());
273 *
for (
unsigned int slice=0; slice<n_slices; ++slice)
275 *
for (
unsigned int i=0; i<input.n_vertices(); ++i)
278 *
if (input.get_used_vertices()[i])
280 *
const Point<2> &v = input.get_vertices()[i];
281 * points[i+slice*input.n_vertices()](0) = v(0);
282 * points[i+slice*input.n_vertices()](1) = v(1);
283 * points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1);
289 * cell = input.begin_active(); cell != input.end(); ++cell)
291 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
294 *
for (
unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
297 * = cell->vertex_index(v)+slice*input.n_used_vertices();
299 * = cell->vertex_index(v)+(slice+1)*input.n_used_vertices();
302 * this_cell.material_id = cell->material_id();
303 * cells.push_back(this_cell);
309 * s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
311 * cell = input.begin_active(); cell != input.end(); ++cell)
314 *
for (
unsigned int f=0; f<4; ++f)
315 *
if (cell->at_boundary(f))
318 * bid =
std::max(bid, quad.boundary_id);
319 *
for (
unsigned int slice=0; slice<n_slices-1; ++slice)
321 * quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_used_vertices();
322 * quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_used_vertices();
323 * quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_used_vertices();
324 * quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_used_vertices();
325 * s.boundary_quads.push_back(quad);
331 * cell = input.begin_active(); cell != input.end(); ++cell)
335 * quad.vertices[0] = cell->vertex_index(0);
336 * quad.vertices[1] = cell->vertex_index(1);
337 * quad.vertices[2] = cell->vertex_index(2);
338 * quad.vertices[3] = cell->vertex_index(3);
339 * s.boundary_quads.push_back(quad);
341 * quad.boundary_id = bid + 2;
342 *
for (
int i=0; i<4; ++i)
343 * quad.vertices[i] += (n_slices-1)*input.n_used_vertices();
344 * s.boundary_quads.push_back(quad);
347 * result.create_triangulation (points,
352 *
namespace Evaluation
364 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
365 * + stress[1][1]*stress[1][1]
366 * - stress[0][0]*stress[1][1]
367 * + 3*stress[0][1]*stress[0][1]);
368 * }
else if (dim == 3)
370 * von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
371 * + stress[1][1]*stress[1][1]
372 * + stress[2][2]*stress[2][2]
373 * - stress[0][0]*stress[1][1]
374 * - stress[1][1]*stress[2][2]
375 * - stress[0][0]*stress[2][2]
376 * + 3*( stress[0][1]*stress[0][1]
377 * +stress[1][2]*stress[1][2]
378 * +stress[0][2]*stress[0][2]) );
383 * -----------------------------------------------
384 *
"Perforated_strip_tension"
386 *
const double von_Mises_stress =
std::sqrt( stress[0][0]*stress[0][0]
387 * + stress[1][1]*stress[1][1]
388 * - stress[0][0]*stress[1][1]
389 * + 3*stress[0][1]*stress[0][1]);
390 * -----------------------------------------------
392 * plane strain / 3
d case
398 * -----------------------------------------------
410 *
return von_Mises_stress;
415 *
class PointValuesEvaluation
418 * PointValuesEvaluation (
const Point<dim> &evaluation_point);
426 * <<
"The evaluation point " << arg1
427 * <<
" was not found among the vertices of the present grid.");
434 * PointValuesEvaluation<dim>::
435 * PointValuesEvaluation (
const Point<dim> &evaluation_point)
437 * evaluation_point (evaluation_point)
444 * PointValuesEvaluation<dim>::
449 *
const unsigned int dofs_per_vertex = dof_handler.get_fe().dofs_per_vertex;
451 * ExcDimensionMismatch (
point_values.size(), dofs_per_vertex));
455 * cell = dof_handler.begin_active(),
456 * endc = dof_handler.end();
457 *
bool evaluation_point_found =
false;
458 *
for (; (cell!=endc) && !evaluation_point_found; ++cell)
460 *
if (cell->is_locally_owned() && !evaluation_point_found)
461 *
for (
unsigned int vertex=0;
462 * vertex<GeometryInfo<dim>::vertices_per_cell;
465 *
if (cell->vertex(vertex).distance (evaluation_point)
467 * cell->diameter() * 1
e-8)
469 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
471 *
point_values[id] = solution(cell->vertex_dof_index(vertex,
id));
474 * evaluation_point_found =
true;
481 * ExcEvaluationPointNotFound(evaluation_point));
490 * <a name=
"elastoplastic.cc-ThecodePointHistorycodeclass"></a>
491 * <h3>The <code>PointHistory</code>
class</h3>
495 * As was mentioned in the introduction, we have to store the old stress in
496 * quadrature
point so that we can compute the residual forces at
this point
497 * during the next time step. This alone would not warrant a structure with
498 * only one member, but in more complicated applications, we would have to
499 * store more information in quadrature points as well, such as the history
500 * variables of plasticity, etc. In essence, we have to store everything
501 * that affects the present state of the material here, which in plasticity
502 * is determined by the deformation history variables.
506 * We will not give
this class any meaningful functionality beyond being
507 * able to store
data, i.e. there are no constructors, destructors, or other
508 * member
functions. In such cases of `dumb
' classes, we usually opt to
509 * declare them as <code>struct</code> rather than <code>class</code>, to
510 * indicate that they are closer to C-style structures than C++-style
515 * struct PointHistory
517 * SymmetricTensor<2,dim> old_stress;
518 * SymmetricTensor<2,dim> old_strain;
526 * <a name="elastoplastic.cc-ThecodeConstitutiveLawcodeclasstemplate"></a>
527 * <h3>The <code>ConstitutiveLaw</code> class template</h3>
531 * This class provides an interface for a constitutive law, i.e., for the
532 * relationship between strain @f$\varepsilon(\mathbf u)@f$ and stress
533 * @f$\sigma@f$. In this example we are using an elastoplastic material behavior
534 * with linear, isotropic hardening. Such materials are characterized by
535 * Young's modulus @f$E@f$, Poisson
's ratio @f$\nu@f$, the initial yield stress
536 * @f$\sigma_0@f$ and the isotropic hardening parameter @f$\gamma@f$. For @f$\gamma =
537 * 0@f$ we obtain perfect elastoplastic behavior.
541 * As explained in the paper that describes this program, the first Newton
542 * steps are solved with a completely elastic material model to avoid having
543 * to deal with both nonlinearities (plasticity and contact) at once. To this
544 * end, this class has a function <code>set_sigma_0()</code> that we use later
545 * on to simply set @f$\sigma_0@f$ to a very large value -- essentially
546 * guaranteeing that the actual stress will not exceed it, and thereby
547 * producing an elastic material. When we are ready to use a plastic model, we
548 * set @f$\sigma_0@f$ back to its proper value, using the same function. As a
549 * result of this approach, we need to leave <code>sigma_0</code> as the only
550 * non-const member variable of this class.
554 * class ConstitutiveLaw
557 * ConstitutiveLaw (const double E,
559 * const double sigma_0,
560 * const double gamma);
563 * set_sigma_0 (double sigma_zero);
566 * get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
567 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
570 * get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
571 * const std::vector<Tensor<2, dim> > &point_hessian,
572 * Tensor<5, dim> &stress_strain_tensor_grad) const;
575 * get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor,
576 * SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
577 * SymmetricTensor<4, dim> &stress_strain_tensor) const;
580 * const double kappa;
583 * const double gamma;
585 * const SymmetricTensor<4, dim> stress_strain_tensor_kappa;
586 * const SymmetricTensor<4, dim> stress_strain_tensor_mu;
591 * The constructor of the ConstitutiveLaw class sets the required material
592 * parameter for our deformable body. Material parameters for elastic
593 * isotropic media can be defined in a variety of ways, such as the pair @f$E,
594 * \nu@f$ (elastic modulus and Poisson's number),
using the Lame parameters
595 * @f$\lambda,mu@f$ or several other commonly used conventions. Here, the
596 * constructor takes a description of material parameters in the form of
597 * @f$E,\nu@f$, but since
this turns out to these are not the coefficients that
598 * appear in the equations of the plastic projector, we immediately convert
599 * them into the more suitable set @f$\kappa,\mu@f$ of bulk and shear moduli. In
600 * addition, the constructor takes @f$\sigma_0@f$ (the yield stress absent any
601 * plastic strain) and @f$\gamma@f$ (the hardening parameter) as arguments. In
602 *
this constructor, we also compute the two principal components of the
603 * stress-strain relation and its linearization.
607 * ConstitutiveLaw<dim>::ConstitutiveLaw (
double E,
614 * --------------------
616 * kappa (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (3 * (1 - 2 * (nu / (1+nu))))),
617 * mu (((E*(1+2*nu)) / (
std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))),
618 * --------------------
619 * 3
d and plane strain
622 * kappa (E / (3 * (1 - 2 * nu))),
623 * mu (E / (2 * (1 + nu))),
626 * --------------------
631 * stress_strain_tensor_kappa (kappa
633 * unit_symmetric_tensor<dim>())),
634 * stress_strain_tensor_mu (2 * mu
635 * * (identity_tensor<dim>()
637 * unit_symmetric_tensor<dim>()) / 3.0))
643 * ConstitutiveLaw<dim>::set_sigma_0 (
double sigma_zero)
645 * sigma_0 = sigma_zero;
652 * <a name=
"elastoplastic.cc-ConstitutiveLawget_stress_strain_tensor"></a>
653 * <h4>ConstitutiveLaw::get_stress_strain_tensor</h4>
657 * This is the principal component of the constitutive law. It projects the
658 * deviatoric part of the stresses in a quadrature
point back to the yield
659 * stress (i.e., the original yield stress @f$\sigma_0@f$ plus the term that
660 * describes linear isotropic hardening). We need
this function to calculate
661 * the nonlinear residual in PlasticityContactProblem::residual_nl_system. The
662 * computations follow the formulas laid out in the introduction.
666 * The function returns whether the quadrature
point is plastic to allow
for
667 * some statistics
downstream on how many of the quadrature points are
668 * plastic and how many are elastic.
673 * ConstitutiveLaw<dim>::
678 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
684 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
687 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
689 * stress_strain_tensor = stress_strain_tensor_mu;
690 *
if (von_Mises_stress > sigma_0)
692 *
const double beta = sigma_0 / von_Mises_stress;
693 * stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
696 * stress_strain_tensor += stress_strain_tensor_kappa;
698 *
return (von_Mises_stress > sigma_0);
704 * ConstitutiveLaw<dim>::
710 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
713 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
715 *
if (von_Mises_stress > sigma_0)
718 *
const double deviator_strain_tensor_norm = deviator_strain_tensor.norm();
719 *
const double multiplier = -(1-
gamma)*sigma_0/(2*mu*
std::pow(deviator_strain_tensor_norm,3));
722 * multiplier_vector = 0;
724 *
for (
unsigned int i=0; i!=dim; ++i)
725 *
for (
unsigned int m=0; m!=dim; ++m)
726 *
for (
unsigned int n=0; n!=dim; ++n)
728 * multiplier_vector(i) += deviator_strain_tensor[m][n] *
729 * ( 0.5*( point_hessian[m][n][i] + point_hessian[n][m][i] )
730 * + ( m==n && dim==2 ? -1/dim*(point_hessian[0][0][i]
731 * + point_hessian[1][1][i]) : 0 )
732 * + ( m==n && dim==3 ? -1/dim*(point_hessian[0][0][i]
733 * + point_hessian[1][1][i]
734 * + point_hessian[2][2][i]) : 0 ) );
739 * -----------------------------------------------
740 *
"Perforated_strip_tension"
743 * -----------------------------------------------
745 * plane strain / 3
d case
748 *
const double VM_factor =
std::sqrt(1.5);
751 * -----------------------------------------------
757 *
for (
unsigned int i=0; i!=dim; ++i)
758 *
for (
unsigned int j=0; j!=dim; ++j)
759 *
for (
unsigned int k=0; k!=dim; ++k)
760 *
for (
unsigned int l=0;
l!=dim; ++
l)
761 *
for (
unsigned int m=0; m!=dim; ++m)
763 * stress_strain_tensor_grad[i][j][k][
l][m] = 1/VM_factor
765 * * stress_strain_tensor_mu[i][j][k][
l]
766 * * multiplier_vector(m);
772 * stress_strain_tensor_grad = 0;
775 *
return (von_Mises_stress > sigma_0);
782 * <a name=
"elastoplastic.cc-ConstitutiveLawget_linearized_stress_strain_tensors"></a>
783 * <h4>ConstitutiveLaw::get_linearized_stress_strain_tensors</h4>
787 * This function returns the linearized stress strain tensor, linearized
788 * around the solution @f$u^{i-1}@f$ of the previous Newton step @f$i-1@f$. The
789 * parameter <code>strain_tensor</code> (commonly denoted
790 * @f$\varepsilon(u^{i-1})@f$) must be passed as an argument, and serves as the
791 * linearization
point. The function returns the derivative of the nonlinear
792 * constitutive law in the variable stress_strain_tensor, as well as the
793 * stress-strain tensor of the linearized problem in
794 * stress_strain_tensor_linearized. See
795 * PlasticityContactProblem::assemble_nl_system where
this function is used.
800 * ConstitutiveLaw<dim>::
806 * stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
809 * stress_strain_tensor = stress_strain_tensor_mu;
810 * stress_strain_tensor_linearized = stress_strain_tensor_mu;
813 *
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
814 *
const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);
816 *
if (von_Mises_stress > sigma_0)
818 *
const double beta = sigma_0 / von_Mises_stress;
819 * stress_strain_tensor *= (
gamma + (1 -
gamma) * beta);
820 * stress_strain_tensor_linearized *= (
gamma + (1 -
gamma) * beta);
821 * deviator_stress_tensor /= deviator_stress_tensor_norm;
822 * stress_strain_tensor_linearized -= (1 -
gamma) * beta * 2 * mu
824 * deviator_stress_tensor);
827 * stress_strain_tensor += stress_strain_tensor_kappa;
828 * stress_strain_tensor_linearized += stress_strain_tensor_kappa;
833 * Finally, below we will need a function that computes the rotation
matrix
834 * induced by a displacement at a given
point. In fact, of course, the
835 * displacement at a single
point only has a direction and a magnitude, it
836 * is the change in direction and magnitude that induces rotations. In
838 * displacement, or, more specifically, from the curl.
842 * The formulas by which the rotation matrices are determined are a little
843 * awkward, especially in 3
d. For 2
d, there is a simpler way, so we
844 * implement
this function twice, once
for 2
d and once
for 3
d, so that we
845 * can compile and use the program in both space dimensions
if so desired --
846 * after all, deal.II is all about dimension
independent programming and
847 * reuse of algorithm thoroughly tested with cheap computations in 2
d,
for
848 * the more expensive computations in 3
d. Here is one
case, where we have to
849 * implement different algorithms
for 2
d and 3
d, but then can write the rest
850 * of the program in a way that is
independent of the space dimension.
854 * So, without further ado to the 2
d implementation:
858 * get_rotation_matrix (
const std::vector<
Tensor<1,2> > &grad_u)
862 * First, compute the curl of the velocity field from the
gradients. Note
863 * that we are in 2
d, so the rotation is a
scalar:
866 *
const double curl = (grad_u[1][0] - grad_u[0][1]);
870 * From
this, compute the
angle of rotation:
877 * And from
this, build the antisymmetric rotation
matrix:
880 *
const double t[2][2] = {{
cos(angle),
sin(angle) },
881 * {-
sin(angle),
cos(angle) }
889 * The 3
d case is a little more contrived:
893 * get_rotation_matrix (
const std::vector<
Tensor<1,3> > &grad_u)
897 * Again
first compute the curl of the velocity field. This time, it is a
901 *
const Point<3> curl (grad_u[2][1] - grad_u[1][2],
902 * grad_u[0][2] - grad_u[2][0],
903 * grad_u[1][0] - grad_u[0][1]);
907 * From
this vector,
using its magnitude, compute the tangent of the
angle
908 * of rotation, and from it the actual
angle:
911 *
const double tan_angle =
std::sqrt(curl*curl);
916 * Now, here
's one problem: if the angle of rotation is too small, that
917 * means that there is no rotation going on (for example a translational
918 * motion). In that case, the rotation matrix is the identity matrix.
922 * The reason why we stress that is that in this case we have that
923 * <code>tan_angle==0</code>. Further down, we need to divide by that
924 * number in the computation of the axis of rotation, and we would get
925 * into trouble when dividing doing so. Therefore, let's shortcut
this and
932 *
static const double rotation[3][3]
933 * = {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } };
940 * Otherwise compute the real rotation
matrix. The algorithm
for this is
941 * not exactly obvious, but can be found in a number of books,
942 * particularly on computer games where rotation is a very frequent
943 * operation. Online, you can find a description at
945 * the signs as here) at
951 *
const double t = 1-c;
953 *
const Point<3> axis = curl/tan_angle;
954 *
const double rotation[3][3]
956 * t *axis[0] *axis[0]+c,
957 * t *axis[0] *axis[1]+s *axis[2],
958 * t *axis[0] *axis[2]-s *axis[1]
961 * t *axis[0] *axis[1]-s *axis[2],
962 * t *axis[1] *axis[1]+c,
963 * t *axis[1] *axis[2]+s *axis[0]
966 * t *axis[0] *axis[2]+s *axis[1],
967 * t *axis[1] *axis[1]-s *axis[0],
968 * t *axis[2] *axis[2]+c
977 * <h3>Equation
data: Body forces, boundary forces,
978 * incremental boundary
values</h3>
982 * The following should be relatively standard. We need classes
for
983 * the boundary forcing term (which we here choose to be zero)
984 * and incremental boundary
values.
987 *
namespace EquationData
1233 * ----------------------------- TimoshenkoBeam ---------------------------------------
1436 * ------------------------- Thick_tube_internal_pressure ----------------------------------
1605 * ------------------------- Perforated_strip_tension ----------------------------------
1775 * ------------------------- Cantiliver_beam_3d ----------------------------------
1778 *
template <
int dim>
1779 *
class IncrementalBoundaryForce :
public Function<dim>
1782 * IncrementalBoundaryForce (
const double present_time,
1783 *
const double end_time);
1795 *
const double present_time,
1801 *
template <
int dim>
1802 * IncrementalBoundaryForce<dim>::
1803 * IncrementalBoundaryForce (
const double present_time,
1804 *
const double end_time)
1807 * present_time (present_time),
1808 * end_time (end_time),
1813 *
template <
int dim>
1815 * IncrementalBoundaryForce<dim>::vector_value (
const Point<dim> &p,
1820 * ExcDimensionMismatch (
values.size(), dim));
1822 *
const double eps = 1.e-7 * height;
1826 * pressure should be imposed on the top surface, y = height
1837 * compute the fraction of imposed force
1840 *
const double frac = present_time/end_time;
1845 *
template <
int dim>
1847 * IncrementalBoundaryForce<dim>::
1848 * vector_value_list (
const std::vector<
Point<dim> > &points,
1851 *
const unsigned int n_points = points.size();
1853 *
Assert (value_list.size() == n_points,
1854 * ExcDimensionMismatch (value_list.size(), n_points));
1856 *
for (
unsigned int p=0; p<n_points; ++p)
1857 * IncrementalBoundaryForce<dim>::vector_value (points[p], value_list[p]);
1861 *
template <
int dim>
1865 * BodyForce () :
Functions::ZeroFunction<dim> (dim) {}
1869 *
template <
int dim>
1870 *
class IncrementalBoundaryValues :
public Function<dim>
1873 * IncrementalBoundaryValues (
const double present_time,
1874 *
const double end_time);
1887 *
const double present_time,
1892 *
template <
int dim>
1893 * IncrementalBoundaryValues<dim>::
1894 * IncrementalBoundaryValues (
const double present_time,
1895 *
const double end_time)
1898 * present_time (present_time),
1899 * end_time (end_time)
1903 *
template <
int dim>
1905 * IncrementalBoundaryValues<dim>::
1910 * ExcDimensionMismatch (
values.size(), dim));
1917 *
template <
int dim>
1919 * IncrementalBoundaryValues<dim>::
1920 * vector_value_list (
const std::vector<
Point<dim> > &points,
1923 *
const unsigned int n_points = points.size();
1925 *
Assert (value_list.size() == n_points,
1926 * ExcDimensionMismatch (value_list.size(), n_points));
1928 *
for (
unsigned int p=0; p<n_points; ++p)
1929 * IncrementalBoundaryValues<dim>::vector_value (points[p], value_list[p]);
1934 * -------------------------------------------------------------------------------
1940 *
namespace DualFunctional
1943 *
template <
int dim>
1951 *
const ConstitutiveLaw<dim> &constitutive_law,
1957 *
template <
int dim>
1958 *
class PointValuesEvaluation :
public DualFunctionalBase<dim>
1961 * PointValuesEvaluation (
const Point<dim> &evaluation_point);
1967 *
const ConstitutiveLaw<dim> &constitutive_law,
1973 * <<
"The evaluation point " << arg1
1974 * <<
" was not found among the vertices of the present grid.");
1981 *
template <
int dim>
1982 * PointValuesEvaluation<dim>::
1983 * PointValuesEvaluation (
const Point<dim> &evaluation_point)
1985 * evaluation_point (evaluation_point)
1989 *
template <
int dim>
1991 * PointValuesEvaluation<dim>::
1994 *
const ConstitutiveLaw<dim> &,
1998 * rhs_dual.reinit (dof_handler_dual.n_dofs());
1999 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2002 * cell_dual = dof_handler_dual.begin_active(),
2003 * endc_dual = dof_handler_dual.end();
2004 *
for (; cell_dual!=endc_dual; ++cell_dual)
2005 *
for (
unsigned int vertex=0;
2006 * vertex<GeometryInfo<dim>::vertices_per_cell;
2008 *
if (cell_dual->vertex(vertex).distance(evaluation_point)
2011 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
2013 * rhs_dual(cell_dual->vertex_dof_index(vertex,
id)) = 1;
2018 *
AssertThrow (
false, ExcEvaluationPointNotFound(evaluation_point));
2022 *
template <
int dim>
2023 *
class PointXDerivativesEvaluation :
public DualFunctionalBase<dim>
2026 * PointXDerivativesEvaluation (
const Point<dim> &evaluation_point);
2032 *
const ConstitutiveLaw<dim> &constitutive_law,
2038 * <<
"The evaluation point " << arg1
2039 * <<
" was not found among the vertices of the present grid.");
2046 *
template <
int dim>
2047 * PointXDerivativesEvaluation<dim>::
2048 * PointXDerivativesEvaluation (
const Point<dim> &evaluation_point)
2050 * evaluation_point (evaluation_point)
2054 *
template <
int dim>
2056 * PointXDerivativesEvaluation<dim>::
2059 *
const ConstitutiveLaw<dim> &,
2063 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2066 *
FEValues<dim> fe_values (dof_handler_dual.get_fe(), quadrature,
2070 *
const unsigned int n_q_points = fe_values.n_quadrature_points;
2071 *
Assert ( n_q_points==quadrature.size() , ExcInternalError() );
2072 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2075 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2077 *
double total_volume = 0;
2080 * cell = dof_handler_dual.begin_active(),
2081 * endc = dof_handler_dual.end();
2082 *
for (; cell!=endc; ++cell)
2083 *
if (cell->center().distance(evaluation_point) <=
2086 * fe_values.reinit (cell);
2089 *
for (
unsigned int q=0; q<n_q_points; ++q)
2091 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2093 * cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
2094 * fe_values.JxW (q);
2097 * total_volume += fe_values.JxW (q);
2100 * cell->get_dof_indices (local_dof_indices);
2101 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2103 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2108 * ExcEvaluationPointNotFound(evaluation_point));
2110 * rhs_dual *= 1./total_volume;
2115 *
template <
int dim>
2116 *
class MeanDisplacementFace :
public DualFunctionalBase<dim>
2119 * MeanDisplacementFace (
const unsigned int face_id,
2120 *
const std::vector<bool> comp_mask);
2126 *
const ConstitutiveLaw<dim> &constitutive_law,
2131 *
const unsigned int face_id;
2132 *
const std::vector<bool> comp_mask;
2136 *
template <
int dim>
2137 * MeanDisplacementFace<dim>::
2138 * MeanDisplacementFace (
const unsigned int face_id,
2139 *
const std::vector<bool> comp_mask )
2141 * face_id (face_id),
2142 * comp_mask (comp_mask)
2145 * ExcDimensionMismatch (comp_mask.size(), dim) );
2149 *
template <
int dim>
2151 * MeanDisplacementFace<dim>::
2154 *
const ConstitutiveLaw<dim> &,
2160 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2162 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2166 *
const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
2167 *
const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
2168 *
const unsigned int n_face_q_points = face_quadrature.size();
2171 * ExcDimensionMismatch (dofs_per_vertex, dim) );
2173 * std::vector<unsigned int> comp_vector(dofs_per_vertex);
2174 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2178 * comp_vector[i] = 1;
2184 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2188 * bound_size :
size of the boundary, in 2
d is the length
2189 * and in the 3
d case, area
2192 *
double bound_size = 0.;
2195 * cell = dof_handler_dual.begin_active(),
2196 * endc = dof_handler_dual.end();
2197 *
bool evaluation_face_found =
false;
2198 *
for (; cell!=endc; ++cell)
2201 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2203 *
if (cell->face(face)->at_boundary()
2205 * cell->face(face)->boundary_id() == face_id)
2207 *
if (!evaluation_face_found)
2209 * evaluation_face_found =
true;
2211 * fe_face_values.reinit (cell, face);
2213 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2215 * bound_size += fe_face_values.JxW(q_point);
2217 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2219 *
const unsigned int
2220 * component_i = dof_handler_dual.get_fe().system_to_component_index(i).first;
2222 * cell_rhs(i) += (fe_face_values.shape_value(i,q_point) *
2223 * comp_vector[component_i] *
2224 * fe_face_values.JxW(q_point));
2232 * cell->get_dof_indices (local_dof_indices);
2233 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2235 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2240 *
AssertThrow(evaluation_face_found, ExcInternalError());
2242 * rhs_dual /= bound_size;
2247 *
template <
int dim>
2248 *
class MeanStressFace :
public DualFunctionalBase<dim>
2251 * MeanStressFace (
const unsigned int face_id,
2252 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2258 *
const ConstitutiveLaw<dim> &constitutive_law,
2263 *
const unsigned int face_id;
2264 *
const std::vector<std::vector<unsigned int> > comp_stress;
2268 *
template <
int dim>
2269 * MeanStressFace<dim>::
2270 * MeanStressFace (
const unsigned int face_id,
2271 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2273 * face_id (face_id),
2274 * comp_stress (comp_stress)
2277 * ExcDimensionMismatch (comp_stress.size(), dim) );
2281 *
template <
int dim>
2283 * MeanStressFace<dim>::
2286 *
const ConstitutiveLaw<dim> &constitutive_law,
2292 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2294 *
const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
2298 *
FEFaceValues<dim> fe_face_values_dual (dof_handler_dual.get_fe(), face_quadrature,
2301 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2302 *
const unsigned int n_face_q_points = face_quadrature.size();
2304 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_face_q_points);
2309 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2313 * bound_size :
size of the boundary, in 2
d is the length
2314 * and in the 3
d case, area
2317 *
double bound_size = 0.;
2319 *
bool evaluation_face_found =
false;
2322 * cell_dual = dof_handler_dual.begin_active(),
2323 * endc_dual = dof_handler_dual.end(),
2324 * cell = dof_handler.begin_active();
2328 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2331 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2333 *
if (cell_dual->face(face)->at_boundary()
2335 * cell_dual->face(face)->boundary_id() == face_id)
2337 *
if (!evaluation_face_found)
2339 * evaluation_face_found =
true;
2342 * fe_face_values.reinit (cell, face);
2343 * fe_face_values_dual.reinit (cell_dual, face);
2345 * fe_face_values[displacement].get_function_symmetric_gradients(solution,
2348 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2350 * bound_size += fe_face_values_dual.JxW(q_point);
2352 * constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2353 * stress_strain_tensor);
2355 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2358 * stress_phi_i = stress_strain_tensor
2359 * * fe_face_values_dual[displacement].symmetric_gradient(i, q_point);
2361 *
for (
unsigned int k=0; k!=dim; ++k)
2363 *
for (
unsigned int l=0;
l!=dim; ++
l)
2365 *
if ( comp_stress[k][l] == 1 )
2367 * cell_rhs(i) += stress_phi_i[k][
l]
2369 * fe_face_values_dual.JxW(q_point);
2382 * cell_dual->get_dof_indices (local_dof_indices);
2383 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2385 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2390 *
AssertThrow(evaluation_face_found, ExcInternalError());
2392 * rhs_dual /= bound_size;
2397 *
template <
int dim>
2398 *
class MeanStressDomain :
public DualFunctionalBase<dim>
2401 * MeanStressDomain (
const std::string &base_mesh,
2402 *
const std::vector<std::vector<unsigned int> > &comp_stress);
2408 *
const ConstitutiveLaw<dim> &constitutive_law,
2413 *
const std::string base_mesh;
2414 *
const std::vector<std::vector<unsigned int> > comp_stress;
2418 *
template <
int dim>
2419 * MeanStressDomain<dim>::
2420 * MeanStressDomain (
const std::string &base_mesh,
2421 *
const std::vector<std::vector<unsigned int> > &comp_stress )
2423 * base_mesh (base_mesh),
2424 * comp_stress (comp_stress)
2427 * ExcDimensionMismatch (comp_stress.size(), dim) );
2431 *
template <
int dim>
2433 * MeanStressDomain<dim>::
2436 *
const ConstitutiveLaw<dim> &constitutive_law,
2440 *
AssertThrow (base_mesh ==
"Cantiliver_beam_3d", ExcNotImplemented());
2445 * Mean stress at the specified domain is of interest.
2446 * The interest domains are located on the bottom and top of the flanges
2447 * close to the clamped face, z = 0
2448 * top domain: height/2 - thickness_flange <= y <= height/2
2449 * 0 <= z <= 2 * thickness_flange
2450 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
2451 * 0 <= z <= 2 * thickness_flange
2457 *
const double height = 200
e-3,
2458 * thickness_flange = 10
e-3;
2460 * rhs_dual.reinit (dof_handler_dual.n_dofs());
2462 *
const QGauss<dim> quadrature_formula(dof_handler_dual.get_fe().tensor_degree()+1);
2464 *
FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
2466 *
FEValues<dim> fe_values_dual (dof_handler_dual.get_fe(), quadrature_formula,
2469 *
const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
2470 *
const unsigned int n_q_points = quadrature_formula.size();
2472 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
2477 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);
2481 * domain_size :
size of the interested domain, in 2
d is the area
2485 *
double domain_size = 0.;
2487 *
bool evaluation_domain_found =
false;
2490 * cell_dual = dof_handler_dual.begin_active(),
2491 * endc_dual = dof_handler_dual.end(),
2492 * cell = dof_handler.begin_active();
2496 *
for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
2498 *
const double y = cell->center()[1],
2499 * z = cell->center()[2];
2502 * top domain: height/2 - thickness_flange <= y <= height/2
2503 * 0 <= z <= 2 * thickness_flange
2504 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
2505 * 0 <= z <= 2 * thickness_flange
2508 *
if ( ((z > 0) && (z < 2*thickness_flange)) &&
2509 * ( ((y > height/2 - thickness_flange) && (y < height/2)) ||
2510 * ((y > -height/2) && (y < -height/2 + thickness_flange)) ) )
2514 *
if (!evaluation_domain_found)
2516 * evaluation_domain_found =
true;
2519 * fe_values.reinit(cell);
2520 * fe_values_dual.reinit(cell_dual);
2522 * fe_values[displacement].get_function_symmetric_gradients(solution,
2525 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
2527 * domain_size += fe_values_dual.JxW(q_point);
2529 * constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
2530 * stress_strain_tensor);
2532 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2535 * stress_phi_i = stress_strain_tensor
2536 * * fe_values_dual[displacement].symmetric_gradient(i, q_point);
2538 *
for (
unsigned int k=0; k!=dim; ++k)
2540 *
for (
unsigned int l=0;
l!=dim; ++
l)
2542 *
if ( comp_stress[k][l] == 1 )
2544 * cell_rhs(i) += stress_phi_i[k][
l]
2546 * fe_values_dual.JxW(q_point);
2558 * cell_dual->get_dof_indices (local_dof_indices);
2559 *
for (
unsigned int i=0; i<dofs_per_cell_dual; ++i)
2561 * rhs_dual(local_dof_indices[i]) += cell_rhs(i);
2566 *
AssertThrow(evaluation_domain_found, ExcInternalError());
2568 * rhs_dual /= domain_size;
2573 *
template <
int dim>
2574 *
class MeanStrainEnergyFace :
public DualFunctionalBase<dim>
2577 * MeanStrainEnergyFace (
const unsigned int face_id,
2581 *
void assemble_rhs_nonlinear (
const DoFHandler<dim> &primal_dof_handler,
2587 *
const unsigned int face_id;
2593 *
template <
int dim>
2594 * MeanStrainEnergyFace<dim>::
2595 * MeanStrainEnergyFace (
const unsigned int face_id,
2599 * face_id (face_id),
2600 * lambda_function (&lambda_function),
2601 * mu_function (&mu_function)
2605 *
template <
int dim>
2607 * MeanStrainEnergyFace<dim>::
2615 * Assemble right hand side of the dual problem when the quantity of interest is
2616 * a nonlinear functional. In
this case, the QoI should be linearized which depends
2617 * on the solution of the primal problem.
2618 * The extractor of the linearized QoI functional is the
gradient of the the original
2619 * QoI functional with the primal solution
values.
2627 * rhs.reinit (dof_handler.n_dofs());
2629 *
const QGauss<dim-1> face_quadrature(dof_handler.get_fe().tensor_degree()+1);
2630 *
FEFaceValues<dim> primal_fe_face_values (primal_dof_handler.get_fe(), face_quadrature,
2638 *
const unsigned int dofs_per_vertex = primal_dof_handler.get_fe().dofs_per_vertex;
2639 *
const unsigned int n_face_q_points = face_quadrature.size();
2640 *
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
2643 * ExcDimensionMismatch (dofs_per_vertex, dim) );
2645 * std::vector< std::vector< Tensor<1,dim> > > primal_solution_gradients;
2646 * primal_solution_gradients.resize(n_face_q_points);
2648 * std::vector<std::vector<Tensor<2,dim> > > primal_solution_hessians;
2649 * primal_solution_hessians.resize (n_face_q_points);
2651 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
2653 * primal_solution_gradients[i].resize (dofs_per_vertex);
2654 * primal_solution_hessians[i].resize (dofs_per_vertex);
2657 * std::vector<double> lambda_values (n_face_q_points);
2658 * std::vector<double> mu_values (n_face_q_points);
2662 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2666 * bound_size :
size of the boundary, in 2
d is the length
2667 * and in the 3
d case, area
2670 *
double bound_size = 0.;
2672 *
bool evaluation_face_found =
false;
2675 * primal_cell = primal_dof_handler.begin_active(),
2676 * primal_endc = primal_dof_handler.end();
2679 * cell = dof_handler.begin_active(),
2680 * endc = dof_handler.end();
2682 *
for (; cell!=endc; ++cell, ++primal_cell)
2685 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
2687 *
if (cell->face(face)->at_boundary()
2689 * cell->face(face)->boundary_id() == face_id)
2691 *
if (!evaluation_face_found)
2693 * evaluation_face_found =
true;
2695 * primal_fe_face_values.reinit (primal_cell, face);
2697 * primal_fe_face_values.get_function_gradients (primal_solution,
2698 * primal_solution_gradients);
2700 * primal_fe_face_values.get_function_hessians (primal_solution,
2701 * primal_solution_hessians);
2703 * lambda_function->value_list (primal_fe_face_values.get_quadrature_points(), lambda_values);
2704 * mu_function->value_list (primal_fe_face_values.get_quadrature_points(), mu_values);
2706 * fe_face_values.reinit (cell, face);
2708 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
2710 * bound_size += primal_fe_face_values.JxW(q_point);
2712 *
for (
unsigned int m=0; m<dofs_per_cell; ++m)
2714 *
const unsigned int
2715 * component_m = dof_handler.get_fe().system_to_component_index(m).first;
2717 *
for (
unsigned int i=0; i!=dofs_per_vertex; ++i)
2719 *
for (
unsigned int j=0; j!=dofs_per_vertex; ++j)
2721 * cell_rhs(m) += fe_face_values.shape_value(m,q_point) *
2723 * lambda_values[q_point] *
2725 * primal_solution_hessians[q_point][i][i][component_m] * primal_solution_gradients[q_point][j][j]
2727 * primal_solution_gradients[q_point][i][i] * primal_solution_hessians[q_point][j][j][component_m]
2730 * mu_values[q_point] *
2732 * 2*primal_solution_hessians[q_point][j][i][component_m] * primal_solution_gradients[q_point][j][i]
2734 * primal_solution_hessians[q_point][i][j][component_m] * primal_solution_gradients[q_point][j][i]
2736 * primal_solution_gradients[q_point][i][j] * primal_solution_hessians[q_point][j][i][component_m]
2739 * primal_fe_face_values.JxW(q_point);
2752 * cell->get_dof_indices (local_dof_indices);
2753 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2755 * rhs(local_dof_indices[i]) += cell_rhs(i);
2760 *
AssertThrow(evaluation_face_found, ExcInternalError());
2762 * rhs *= 1./(2*bound_size);
2775 *
template <
int dim>
2782 *
const ConstitutiveLaw<dim> &constitutive_law,
2783 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2784 *
const unsigned int ×tep_no,
2785 *
const std::string &output_dir,
2786 *
const std::string &base_mesh,
2787 *
const double &present_time,
2788 *
const double &end_time);
2790 *
void compute_error_DWR (
Vector<float> &estimated_error_per_cell);
2795 *
void setup_system ();
2796 *
void compute_dirichlet_constraints ();
2797 *
void assemble_matrix ();
2798 *
void assemble_rhs ();
2800 *
void output_results ();
2806 *
const unsigned int fe_degree;
2809 *
const unsigned int fe_degree_dual;
2814 *
const QGauss<dim - 1> face_quadrature_formula;
2824 *
const ConstitutiveLaw<dim> constitutive_law;
2829 *
unsigned int timestep_no;
2830 * std::string output_dir;
2831 *
const std::string base_mesh;
2832 *
double present_time;
2842 *
const ConstitutiveLaw<dim> &constitutive_law,
2843 *
const DualFunctional::DualFunctionalBase<dim> &dual_functional,
2844 *
const unsigned int ×tep_no,
2845 *
const std::string &output_dir,
2846 *
const std::string &base_mesh,
2847 *
const double &present_time,
2848 *
const double &end_time)
2852 * solution(solution),
2853 * fe_degree(fe.tensor_degree()),
2854 * fe_degree_dual(fe_degree + 1),
2855 * fe_dual(
FE_Q<dim>(fe_degree_dual), dim),
2857 * quadrature_formula (fe_degree_dual + 1),
2858 * face_quadrature_formula (fe_degree_dual + 1),
2859 * constitutive_law (constitutive_law),
2861 * dual_functional (&dual_functional),
2862 * timestep_no (timestep_no),
2863 * output_dir (output_dir),
2864 * base_mesh (base_mesh),
2865 * present_time (present_time),
2866 * end_time (end_time)
2871 * DualSolver<dim>::~DualSolver()
2873 * dof_handler_dual.
clear ();
2878 *
void DualSolver<dim>::setup_system()
2880 * dof_handler.distribute_dofs(fe);
2882 * dof_handler_dual.distribute_dofs (fe_dual);
2883 * std::cout <<
" Number of degrees of freedom in dual problem: "
2884 * << dof_handler_dual.n_dofs()
2887 * constraints_hanging_nodes_dual.clear ();
2889 * constraints_hanging_nodes_dual);
2890 * constraints_hanging_nodes_dual.close ();
2892 * compute_dirichlet_constraints();
2894 * sparsity_pattern_dual.reinit (dof_handler_dual.n_dofs(),
2895 * dof_handler_dual.n_dofs(),
2896 * dof_handler_dual.max_couplings_between_dofs());
2901 * constraints_hanging_nodes_dual.condense (sparsity_pattern_dual);
2904 * constraints_dirichlet_and_hanging_nodes_dual.condense (sparsity_pattern_dual);
2906 * sparsity_pattern_dual.compress();
2908 * system_matrix_dual.reinit (sparsity_pattern_dual);
2910 * solution_dual.reinit (dof_handler_dual.n_dofs());
2911 * system_rhs_dual.reinit (dof_handler_dual.n_dofs());
2916 *
void DualSolver<dim>::compute_dirichlet_constraints()
2918 * constraints_dirichlet_and_hanging_nodes_dual.clear ();
2919 * constraints_dirichlet_and_hanging_nodes_dual.merge(constraints_hanging_nodes_dual);
2921 * std::vector<bool> component_mask(dim);
2923 *
if (base_mesh ==
"Timoshenko beam")
2927 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2928 * constraints_dirichlet_and_hanging_nodes_dual,
2931 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
2935 * the boundary x = 0
2938 * component_mask[0] =
true;
2939 * component_mask[1] =
false;
2942 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2943 * constraints_dirichlet_and_hanging_nodes_dual,
2947 * the boundary y = 0
2950 * component_mask[0] =
false;
2951 * component_mask[1] =
true;
2954 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2955 * constraints_dirichlet_and_hanging_nodes_dual,
2958 *
else if (base_mesh ==
"Perforated_strip_tension")
2962 * the boundary x = 0
2965 * component_mask[0] =
true;
2966 * component_mask[1] =
false;
2967 * component_mask[2] =
false;
2970 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2971 * constraints_dirichlet_and_hanging_nodes_dual,
2975 * the boundary y = 0
2978 * component_mask[0] =
false;
2979 * component_mask[1] =
true;
2980 * component_mask[2] =
false;
2983 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2984 * constraints_dirichlet_and_hanging_nodes_dual,
2988 * the boundary y = imposed incremental displacement
2991 * component_mask[0] =
false;
2992 * component_mask[1] =
true;
2993 * component_mask[2] =
false;
2996 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
2997 * constraints_dirichlet_and_hanging_nodes_dual,
3000 *
else if (base_mesh ==
"Cantiliver_beam_3d")
3004 * the boundary x = y = z = 0
3007 * component_mask[0] =
true;
3008 * component_mask[1] =
true;
3009 * component_mask[2] =
true;
3012 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
3013 * constraints_dirichlet_and_hanging_nodes_dual,
3021 * constraints_dirichlet_and_hanging_nodes_dual.close();
3026 *
void DualSolver<dim>::assemble_matrix()
3033 *
const unsigned int dofs_per_cell_dual = fe_dual.dofs_per_cell;
3034 *
const unsigned int n_q_points = quadrature_formula.size();
3038 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell_dual);
3041 * cell_dual = dof_handler_dual.begin_active(),
3042 * endc_dual = dof_handler_dual.end(),
3043 * cell = dof_handler.begin_active();
3047 *
for (; cell_dual != endc_dual; ++cell_dual, ++cell)
3048 *
if (cell_dual->is_locally_owned())
3050 * fe_values.reinit(cell);
3052 * fe_values_dual.reinit(cell_dual);
3055 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3056 * fe_values[displacement].get_function_symmetric_gradients(solution,
3059 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3063 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3064 * stress_strain_tensor_linearized,
3065 * stress_strain_tensor);
3067 *
for (
unsigned int i = 0; i < dofs_per_cell_dual; ++i)
3070 * stress_phi_i = stress_strain_tensor_linearized
3071 * * fe_values_dual[displacement].symmetric_gradient(i, q_point);
3073 *
for (
unsigned int j = 0; j < dofs_per_cell_dual; ++j)
3075 * * fe_values_dual[displacement].symmetric_gradient(j, q_point)
3076 * * fe_values_dual.JxW(q_point));
3082 * cell_dual->get_dof_indices(local_dof_indices);
3083 * constraints_dirichlet_and_hanging_nodes_dual.distribute_local_to_global(cell_matrix,
3084 * local_dof_indices,
3085 * system_matrix_dual);
3093 *
void DualSolver<dim>::assemble_rhs()
3095 * dual_functional->assemble_rhs (dof_handler, solution, constitutive_law,
3096 * dof_handler_dual, system_rhs_dual);
3097 * constraints_dirichlet_and_hanging_nodes_dual.condense (system_rhs_dual);
3102 *
void DualSolver<dim>::solve()
3106 * +++ direct solver +++++++++
3114 * After the decomposition, we can use A_direct like a
matrix representing
3115 * the inverse of our system
matrix, so to compute the solution we just
3116 * have to multiply with the right hand side vector:
3119 * A_direct.vmult(solution_dual, system_rhs_dual);
3123 * ++++ iterative solver ++ CG ++++ doesn
't work
3124 * SolverControl solver_control (5000, 1e-12);
3125 * SolverCG<> cg (solver_control);
3129 * PreconditionSSOR<> preconditioner;
3130 * preconditioner.initialize(system_matrix_dual, 1.2);
3134 * cg.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3139 * ++++ iterative solver ++ BiCGStab ++++++ doesn't work
3146 * preconditioner.
initialize(system_matrix_dual, 1.0);
3150 * bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual,
3155 * +++++++++++++++++++++++++++++++++++++++++++++++++
3161 * constraints_dirichlet_and_hanging_nodes_dual.distribute (solution_dual);
3165 *
void DualSolver<dim>::output_results()
3167 * std::string filename = (output_dir +
"dual-solution-" +
3169 * std::ofstream output (filename.c_str());
3172 * std::vector<std::string> solution_names;
3176 * solution_names.push_back (
"displacement");
3179 * solution_names.push_back (
"x_displacement");
3180 * solution_names.push_back (
"y_displacement");
3183 * solution_names.push_back (
"x_displacement");
3184 * solution_names.push_back (
"y_displacement");
3185 * solution_names.push_back (
"z_displacement");
3188 *
Assert (
false, ExcNotImplemented());
3190 * data_out.add_data_vector (solution_dual, solution_names);
3191 * data_out.build_patches ();
3192 * data_out.write_vtk (output);
3196 *
void DualSolver<dim>::compute_error_DWR (
Vector<float> &estimated_error_per_cell)
3199 * ExcDimensionMismatch (estimated_error_per_cell.size(),
triangulation->n_global_active_cells()));
3203 * solve the dual problem
3207 * assemble_matrix ();
3210 * output_results ();
3214 * compuate the dual weights
3221 * constraints_dirichlet_and_hanging_nodes_dual,
3226 * constraints_hanging_nodes);
3227 * constraints_hanging_nodes.close();
3230 * constraints_dirichlet_and_hanging_nodes_dual,
3233 * constraints_hanging_nodes,
3238 * estimate the error
3248 *
const unsigned int n_q_points = quadrature_formula.size();
3249 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
3253 * std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
3254 *
for (
unsigned int i=0; i!=n_q_points; ++i)
3256 * cell_hessians[i].resize (dim);
3258 * std::vector<Vector<double> > dual_weights_cell_values (n_q_points,
Vector<double>(dim));
3260 *
const EquationData::BodyForce<dim> body_force;
3261 * std::vector<Vector<double> > body_force_values (n_q_points,
Vector<double>(dim));
3271 * fe_face_values_neighbor (fe_dual, face_quadrature_formula,
3279 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
3280 * std::vector<Vector<double> > jump_residual (n_face_q_points,
Vector<double>(dim));
3281 * std::vector<Vector<double> > dual_weights_face_values (n_face_q_points,
Vector<double>(dim));
3283 * std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
3284 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3286 * cell_grads[i].resize (dim);
3288 * std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
3289 *
for (
unsigned int i=0; i!=n_face_q_points; ++i)
3291 * neighbor_grads[i].resize (dim);
3299 *
typename std::map<typename DoFHandler<dim>::face_iterator,
Vector<double> >
3302 * cell = dof_handler_dual.begin_active(),
3303 * endc = dof_handler_dual.
end();
3304 *
for (; cell!=endc; ++cell)
3305 *
if (cell->is_locally_owned())
3307 *
for (
unsigned int face_no=0;
3308 * face_no<GeometryInfo<dim>::faces_per_cell;
3311 * face_integrals[cell->face(face_no)].reinit (dim);
3312 * face_integrals[cell->face(face_no)] = -1e20;
3316 * std::vector<Vector<float> > error_indicators_vector;
3317 * error_indicators_vector.resize(
triangulation->n_active_cells(),
3322 * ----------------- estimate_some -------------------------
3325 * cell = dof_handler_dual.begin_active();
3326 *
unsigned int present_cell = 0;
3327 *
for (; cell!=endc; ++cell, ++present_cell)
3328 *
if (cell->is_locally_owned())
3332 * --------------- integrate_over_cell -------------------
3335 * fe_values.reinit(cell);
3336 * body_force.vector_value_list(fe_values.get_quadrature_points(),
3337 * body_force_values);
3338 * fe_values[displacement].get_function_symmetric_gradients(primal_solution,
3340 * fe_values.get_function_hessians(primal_solution, cell_hessians);
3342 * fe_values.get_function_values(dual_weights,
3343 * dual_weights_cell_values);
3345 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
3347 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
3348 * stress_strain_tensor_linearized,
3349 * stress_strain_tensor);
3350 * constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
3351 * cell_hessians[q_point],
3352 * stress_strain_tensor_grad);
3354 *
for (
unsigned int i=0; i!=dim; ++i)
3356 * error_indicators_vector[present_cell](i) +=
3357 * body_force_values[q_point](i)*
3358 * dual_weights_cell_values[q_point](i)*
3359 * fe_values.JxW(q_point);
3360 *
for (
unsigned int j=0; j!=dim; ++j)
3362 *
for (
unsigned int k=0; k!=dim; ++k)
3364 *
for (
unsigned int l=0;
l!=dim; ++
l)
3366 * error_indicators_vector[present_cell](i) +=
3367 * ( stress_strain_tensor[i][j][k][l]*
3368 * 0.5*(cell_hessians[q_point][k][l][j]
3370 * cell_hessians[q_point][l][k][j])
3371 * + stress_strain_tensor_grad[i][j][k][
l][j] * strain_tensor[q_point][k][
l]
3373 * dual_weights_cell_values[q_point](i) *
3374 * fe_values.JxW(q_point);
3384 * -------------------------------------------------------
3385 * compute face_integrals
3388 *
for (
unsigned int face_no=0;
3389 * face_no<GeometryInfo<dim>::faces_per_cell;
3392 *
if (cell->face(face_no)->at_boundary())
3394 *
for (
unsigned int id=0;
id!=dim; ++id)
3396 * face_integrals[cell->face(face_no)](id) = 0;
3401 *
if ((cell->neighbor(face_no)->has_children() ==
false) &&
3402 * (cell->neighbor(face_no)->level() == cell->level()) &&
3403 * (cell->neighbor(face_no)->index() < cell->index()))
3406 *
if (cell->at_boundary(face_no) ==
false)
3407 *
if (cell->neighbor(face_no)->level() < cell->level())
3411 *
if (cell->face(face_no)->has_children() ==
false)
3415 * ------------- integrate_over_regular_face -----------
3418 * fe_face_values_cell.reinit(cell, face_no);
3419 * fe_face_values_cell.get_function_gradients (primal_solution,
3423 * ExcInternalError());
3424 *
const unsigned int
3425 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3427 * neighbor = cell->neighbor(face_no);
3429 * fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
3430 * fe_face_values_neighbor.get_function_gradients (primal_solution,
3433 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3435 * q_cell_strain_tensor = 0.;
3436 * q_neighbor_strain_tensor = 0.;
3437 *
for (
unsigned int i=0; i!=dim; ++i)
3439 *
for (
unsigned int j=0; j!=dim; ++j)
3441 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3442 * cell_grads[q_point][j][i] );
3443 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3444 * neighbor_grads[q_point][j][i] );
3448 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3449 * cell_stress_strain_tensor);
3450 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3451 * neighbor_stress_strain_tensor);
3453 * jump_residual[q_point] = 0.;
3454 *
for (
unsigned int i=0; i!=dim; ++i)
3456 *
for (
unsigned int j=0; j!=dim; ++j)
3458 *
for (
unsigned int k=0; k!=dim; ++k)
3460 *
for (
unsigned int l=0;
l!=dim; ++
l)
3462 * jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
3463 * q_cell_strain_tensor[k][l]
3465 * neighbor_stress_strain_tensor[i][j][k][l]*
3466 * q_neighbor_strain_tensor[k][l] )*
3467 * fe_face_values_cell.normal_vector(q_point)[j];
3475 * fe_face_values_cell.get_function_values (dual_weights,
3476 * dual_weights_face_values);
3479 * face_integral_vector = 0;
3480 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3482 *
for (
unsigned int i=0; i!=dim; ++i)
3484 * face_integral_vector(i) += jump_residual[q_point](i) *
3485 * dual_weights_face_values[q_point](i) *
3486 * fe_face_values_cell.JxW(q_point);
3490 *
Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
3491 * ExcInternalError());
3493 *
for (
unsigned int i=0; i!=dim; ++i)
3495 *
Assert (face_integrals[cell->face(face_no)](i) == -1e20,
3496 * ExcInternalError());
3497 * face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
3503 * -----------------------------------------------------
3511 * ------------- integrate_over_irregular_face ---------
3515 * face = cell->face(face_no);
3517 * neighbor = cell->neighbor(face_no);
3519 * ExcInternalError());
3520 *
Assert (neighbor->has_children(),
3521 * ExcInternalError());
3523 *
const unsigned int
3524 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
3526 *
for (
unsigned int subface_no=0;
3527 * subface_no<face->n_children(); ++subface_no)
3530 * neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
3531 *
Assert (neighbor_child->face(neighbor_neighbor) ==
3532 * cell->face(face_no)->child(subface_no),
3533 * ExcInternalError());
3535 * fe_subface_values_cell.reinit (cell, face_no, subface_no);
3536 * fe_subface_values_cell.get_function_gradients (primal_solution,
3538 * fe_face_values_neighbor.reinit (neighbor_child,
3539 * neighbor_neighbor);
3540 * fe_face_values_neighbor.get_function_gradients (primal_solution,
3543 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3545 * q_cell_strain_tensor = 0.;
3546 * q_neighbor_strain_tensor = 0.;
3547 *
for (
unsigned int i=0; i!=dim; ++i)
3549 *
for (
unsigned int j=0; j!=dim; ++j)
3551 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
3552 * cell_grads[q_point][j][i] );
3553 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
3554 * neighbor_grads[q_point][j][i] );
3558 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
3559 * cell_stress_strain_tensor);
3560 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
3561 * neighbor_stress_strain_tensor);
3563 * jump_residual[q_point] = 0.;
3564 *
for (
unsigned int i=0; i!=dim; ++i)
3566 *
for (
unsigned int j=0; j!=dim; ++j)
3568 *
for (
unsigned int k=0; k!=dim; ++k)
3570 *
for (
unsigned int l=0;
l!=dim; ++
l)
3572 * jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
3573 * q_cell_strain_tensor[k][l]
3575 * neighbor_stress_strain_tensor[i][j][k][l]*
3576 * q_neighbor_strain_tensor[k][l] )*
3577 * fe_face_values_neighbor.normal_vector(q_point)[j];
3585 * fe_face_values_neighbor.get_function_values (dual_weights,
3586 * dual_weights_face_values);
3589 * face_integral_vector = 0;
3590 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
3592 *
for (
unsigned int i=0; i!=dim; ++i)
3594 * face_integral_vector(i) += jump_residual[q_point](i) *
3595 * dual_weights_face_values[q_point](i) *
3596 * fe_face_values_neighbor.JxW(q_point);
3600 *
for (
unsigned int i=0; i!=dim; ++i)
3602 * face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
3609 *
for (
unsigned int subface_no=0;
3610 * subface_no<face->n_children(); ++subface_no)
3612 *
Assert (face_integrals.find(face->child(subface_no)) !=
3613 * face_integrals.end(),
3614 * ExcInternalError());
3615 *
for (
unsigned int i=0; i!=dim; ++i)
3617 *
Assert (face_integrals[face->child(subface_no)](i) != -1e20,
3618 * ExcInternalError());
3619 *
sum(i) += face_integrals[face->child(subface_no)](i);
3622 *
for (
unsigned int i=0; i!=dim; ++i)
3624 * face_integrals[face](i) =
sum(i);
3630 * -----------------------------------------------------
3640 * ----------------------------------------------------------
3647 * cell = dof_handler_dual.begin_active();
3648 *
for (; cell!=endc; ++cell, ++present_cell)
3649 *
if (cell->is_locally_owned())
3651 *
for (
unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
3654 *
Assert(face_integrals.find(cell->face(face_no)) !=
3655 * face_integrals.end(),
3656 * ExcInternalError());
3658 *
for (
unsigned int id=0;
id!=dim; ++id)
3660 * error_indicators_vector[present_cell](id)
3661 * -= 0.5*face_integrals[cell->face(face_no)](id);
3666 * estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
3676 * <a name=
"elastoplastic.cc-ThecodePlasticityContactProblemcodeclasstemplate"></a>
3677 * <h3>The <code>PlasticityContactProblem</code>
class template</h3>
3681 * This is the main
class of this program and supplies all
functions
3682 * and variables needed to describe
3683 * the nonlinear contact problem. It is
3684 * close to @ref step_41
"step-41" but with some additional
3685 * features like handling hanging nodes,
3686 * a Newton method,
using Trilinos and p4est
3687 *
for parallel distributed computing.
3688 * To deal with hanging nodes makes
3689 * life a bit more complicated since
3691 * We create a Newton method
for the
3692 * active set method
for the contact
3693 * situation and to handle the nonlinear
3694 *
operator for the constitutive law.
3698 * The
general layout of
this class is very much like for most other tutorial programs.
3699 * To make our life a bit easier,
this class reads a set of input parameters from an input file. These
3700 * parameters,
using the
ParameterHandler class, are declared in the <code>declare_parameters</code>
3701 * function (which is
static so that it can be called before we even create an
object of the current
3702 * type), and a
ParameterHandler object that has been used to read an input file will then be passed
3703 * to the constructor of
this class.
3707 * The remaining member
functions are by and large as we have seen in several of the other tutorial
3708 * programs, though with additions
for the current nonlinear system. We will comment on their purpose
3709 * as we get to them further below.
3712 *
template <
int dim>
3713 *
class ElastoPlasticProblem
3723 *
void make_grid ();
3724 *
void setup_system ();
3725 *
void compute_dirichlet_constraints ();
3729 *
void solve_newton_system ();
3730 *
void solve_newton ();
3731 *
void compute_error ();
3733 *
void refine_grid ();
3735 *
void output_results (
const std::string &filename_base);
3739 * Next are three
functions that handle the history variables stored in each
3740 * quadrature
point. The
first one is called before the
first timestep to
3741 * set up a pristine state
for the history variables. It only works on
3742 * those quadrature points on cells that belong to the present processor:
3745 *
void setup_quadrature_point_history ();
3749 * The
second one updates the history variables at the
end of each
3753 *
void update_quadrature_point_history ();
3757 * As far as member variables are concerned, we start with ones that we use to
3758 * indicate the
MPI universe
this program runs on, and then two
numbers
3759 * telling us how many participating processors there are, and where in
3760 *
this world we are., a stream we use to let
3761 * exactly one processor produce output to the console (see @ref step_17
"step-17") and
3762 * a variable that is used to time the various sections of the program:
3773 * The next
group describes the mesh and the finite element space.
3774 * In particular,
for this parallel program, the finite element
3775 * space has associated with it variables that indicate which degrees
3776 * of freedom live on the current processor (the index sets, see
3777 * also @ref step_40
"step-40" and the @ref distributed documentation module) as
3778 * well as a variety of constraints: those imposed by hanging nodes,
3779 * by Dirichlet boundary conditions, and by the active set of
3781 * here, the
first only contains hanging node constraints, the
3782 *
second also those associated with Dirichlet boundary conditions,
3783 * and the third these plus the contact constraints.
3787 * The variable <code>active_set</code> consists of those degrees
3788 * of freedom constrained by the contact, and we use
3789 * <code>fraction_of_plastic_q_points_per_cell</code> to keep
3790 * track of the fraction of quadrature points on each cell where
3791 * the stress equals the yield stress. The latter is only used to
3792 * create graphical output showing the plastic zone, but not
for
3793 * any further computation; the variable is a member variable of
3794 *
this class since the information is computed as a by-product
3795 * of computing the residual, but is used only much later. (Note
3796 * that the vector is a vector of length
equal to the number of
3797 * active cells on the <i>local mesh</i>; it is never used to
3798 * exchange information between processors and can therefore be
3799 * a regular deal.II vector.)
3802 *
const unsigned int n_initial_global_refinements;
3805 *
const unsigned int fe_degree;
3819 * One difference of
this program is that we declare the quadrature
3820 * formula in the
class declaration. The reason is that in all the other
3821 * programs, it didn
't do much harm if we had used different quadrature
3822 * formulas when computing the matrix and the right hand side, for
3823 * example. However, in the present case it does: we store information in
3824 * the quadrature points, so we have to make sure all parts of the program
3825 * agree on where they are and how many there are on each cell. Thus, let
3826 * us first declare the quadrature formula that will be used throughout...
3829 * const QGauss<dim> quadrature_formula;
3830 * const QGauss<dim - 1> face_quadrature_formula;
3834 * ... and then also have a vector of history objects, one per quadrature
3835 * point on those cells for which we are responsible (i.e. we don't store
3836 * history
data for quadrature points on cells that are owned by other
3840 * std::vector<PointHistory<dim> > quadrature_point_history;
3844 * The way
this object is accessed is through a <code>user pointer</code>
3845 * that each cell, face, or edge holds: it is a <code>
void*</code> pointer
3846 * that can be used by application programs to associate arbitrary
data to
3847 * cells, faces, or edges. What the program actually does with
this data
3848 * is within its own responsibility, the library just allocates some space
3849 *
for these pointers, and application programs can set and read the
3850 * pointers
for each of these objects.
3857 * The next block of variables corresponds to the solution
3858 * and the linear systems we need to form. In particular,
this
3859 * includes the Newton
matrix and right hand side; the vector
3860 * that corresponds to the residual (i.e., the Newton right hand
3861 * side) but from which we have not eliminated the various
3862 * constraints and that is used to determine which degrees of
3863 * freedom need to be constrained in the next iteration; and
3864 * a vector that corresponds to the
diagonal of the @f$B@f$
matrix
3865 * briefly mentioned in the introduction and discussed in the
3866 * accompanying paper.
3878 * The next block of variables is then related to the time dependent
3879 * nature of the problem: they denote the length of the time interval
3880 * which we want to simulate, the present time and number of time step,
3881 * and length of present timestep:
3884 *
double present_time;
3885 *
double present_timestep;
3887 *
unsigned int timestep_no;
3891 * The next block contains the variables that describe the material
3895 *
const double e_modulus, nu, sigma_0,
gamma;
3896 * ConstitutiveLaw<dim> constitutive_law;
3900 * And then there is an assortment of other variables that are used
3901 * to identify the mesh we are asked to build as selected by the
3902 * parameter file, the obstacle that is being pushed into the
3903 * deformable body, the mesh refinement strategy, whether to transfer
3904 * the solution from one mesh to the next, and how many mesh
3905 * refinement cycles to perform. As possible, we mark these kinds
3906 * of variables as <code>
const</code> to help the reader identify
3907 * which ones may or may not be modified later on (the output directory
3908 * being an exception -- it is never modified outside the constructor
3909 * but it is awkward to initialize in the member-initializer-list
3910 * following the colon in the constructor since there we have only
3911 * one shot at setting it; the same is
true for the mesh refinement
3915 * const
std::string base_mesh;
3917 *
struct RefinementStrategy
3922 * refine_percentage,
3926 *
typename RefinementStrategy::value refinement_strategy;
3928 *
struct ErrorEstimationStrategy
3934 * weighted_residual_error,
3935 * weighted_kelly_error
3938 *
typename ErrorEstimationStrategy::value error_estimation_strategy;
3942 *
const bool transfer_solution;
3943 * std::string output_dir;
3948 *
unsigned int current_refinement_cycle;
3950 *
const double max_relative_error;
3951 *
float relative_error;
3953 *
const bool show_stresses;
3960 * <a name=
"elastoplastic.cc-ImplementationofthecodePlasticityContactProblemcodeclass"></a>
3961 * <h3>Implementation of the <code>PlasticityContactProblem</code>
class</h3>
3966 * <a name=
"elastoplastic.cc-PlasticityContactProblemdeclare_parameters"></a>
3967 * <h4>PlasticityContactProblem::declare_parameters</h4>
3971 * Let us start with the declaration of
run-time parameters that can be
3972 * selected in the input file. These
values will be read back in the
3973 * constructor of
this class to initialize the member variables of this
3977 *
template <
int dim>
3981 * prm.declare_entry(
"polynomial degree",
"1",
3983 *
"Polynomial degree of the FE_Q finite element space, typically 1 or 2.");
3984 * prm.declare_entry(
"number of initial refinements",
"2",
3986 *
"Number of initial global mesh refinement steps before "
3987 *
"the first computation.");
3988 * prm.declare_entry(
"refinement strategy",
"percentage",
3990 *
"Mesh refinement strategy:\n"
3991 *
" global: one global refinement\n"
3992 *
" percentage: a fixed percentage of cells gets refined using the selected error estimator.");
3993 * prm.declare_entry(
"error estimation strategy",
"kelly_error",
3995 *
"Error estimation strategy:\n"
3996 *
" kelly_error: Kelly error estimator\n"
3997 *
" residual_error: residual-based error estimator\n"
3998 *
" weighted_residual_error: dual weighted residual (Goal-oriented) error estimator.\n");
3999 * prm.declare_entry(
"maximum relative error",
"0.05",
4001 *
"maximum relative error which plays the role of a criteria for refinement.");
4002 * prm.declare_entry(
"number of cycles",
"5",
4004 *
"Number of adaptive mesh refinement cycles to run.");
4005 * prm.declare_entry(
"output directory",
"",
4007 *
"Directory for output files (graphical output and benchmark "
4008 *
"statistics). If empty, use the current directory.");
4009 * prm.declare_entry(
"transfer solution",
"true",
4011 *
"Whether the solution should be used as a starting guess "
4012 *
"for the next finer mesh. If false, then the iteration starts at "
4013 *
"zero on every mesh.");
4014 * prm.declare_entry(
"base mesh",
"Thick_tube_internal_pressure",
4016 *
"Perforated_strip_tension|Cantiliver_beam_3d"),
4017 *
"Select the shape of the domain: 'box' or 'half sphere'");
4018 * prm.declare_entry(
"elasticity modulus",
"2.e11",
4020 *
"Elasticity modulus of the material in MPa (N/mm2)");
4021 * prm.declare_entry(
"Poissons ratio",
"0.3",
4023 *
"Poisson's ratio of the material");
4024 * prm.declare_entry(
"yield stress",
"2.e11",
4026 *
"Yield stress of the material in MPa (N/mm2)");
4027 * prm.declare_entry(
"isotropic hardening parameter",
"0.",
4029 *
"Isotropic hardening parameter of the material");
4030 * prm.declare_entry(
"show stresses",
"false",
4032 *
"Whether illustrates the stresses and von Mises stresses or not.");
4041 * <a name=
"elastoplastic.cc-ThecodePlasticityContactProblemcodeconstructor"></a>
4042 * <h4>The <code>PlasticityContactProblem</code> constructor</h4>
4046 * Given the declarations of member variables as well as the
4047 * declarations of
run-time parameters that are read from the input
4048 * file, there is
nothing surprising in
this constructor. In the body
4049 * we initialize the mesh refinement strategy and the output directory,
4050 * creating such a directory
if necessary.
4053 *
template <
int dim>
4054 * ElastoPlasticProblem<dim>::
4057 * mpi_communicator(MPI_COMM_WORLD),
4060 * pcout(std::cout, this_mpi_process == 0),
4064 * n_initial_global_refinements (prm.get_integer(
"number of initial refinements")),
4066 * fe_degree (prm.get_integer(
"polynomial degree")),
4069 * quadrature_formula (fe_degree + 1),
4070 * face_quadrature_formula (fe_degree + 1),
4072 * e_modulus (prm.get_double(
"elasticity modulus")),
4073 * nu (prm.get_double(
"Poissons ratio")),
4074 * sigma_0(prm.get_double(
"yield stress")),
4075 *
gamma (prm.get_double(
"isotropic hardening parameter")),
4076 * constitutive_law (e_modulus,
4081 * base_mesh (prm.get(
"base mesh")),
4083 * transfer_solution (prm.get_bool(
"transfer solution")),
4085 * table_results_2(),
4086 * table_results_3(),
4087 * max_relative_error (prm.get_double(
"maximum relative error")),
4088 * show_stresses (prm.get_bool(
"show stresses"))
4090 * std::string strat = prm.get(
"refinement strategy");
4091 *
if (strat ==
"global")
4092 * refinement_strategy = RefinementStrategy::refine_global;
4093 *
else if (strat ==
"percentage")
4094 * refinement_strategy = RefinementStrategy::refine_percentage;
4098 * strat = prm.get(
"error estimation strategy");
4099 *
if (strat ==
"kelly_error")
4100 * error_estimation_strategy = ErrorEstimationStrategy::kelly_error;
4101 *
else if (strat ==
"residual_error")
4102 * error_estimation_strategy = ErrorEstimationStrategy::residual_error;
4103 *
else if (strat ==
"weighted_residual_error")
4104 * error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error;
4108 * output_dir = prm.get(
"output directory");
4109 *
if (output_dir !=
"" && *(output_dir.rbegin()) !=
'/')
4110 * output_dir +=
"/";
4111 * mkdir(output_dir.c_str(), 0777);
4113 * pcout <<
" Using output directory '" << output_dir <<
"'" << std::endl;
4114 * pcout <<
" FE degree " << fe_degree << std::endl;
4115 * pcout <<
" transfer solution "
4116 * << (transfer_solution ?
"true" :
"false") << std::endl;
4124 * <a name=
"elastoplastic.cc-PlasticityContactProblemmake_grid"></a>
4125 * <h4>PlasticityContactProblem::make_grid</h4>
4129 * The next block deals with constructing the starting mesh.
4130 * We will use the following helper function and the
first
4131 * block of the <code>make_grid()</code> to construct a
4132 * mesh that corresponds to a half sphere. deal.II has a function
4133 * that creates such a mesh, but it is in the wrong location
4134 * and facing the wrong direction, so we need to
shift and
rotate
4135 * it a bit before
using it.
4139 * For later reference, as described in the documentation of
4141 * has boundary indicator zero,
while the remainder has boundary
4146 * rotate_half_sphere (
const Point<3> &in)
4148 *
return Point<3>(in(2), in(1), -in(0));
4151 *
template <
int dim>
4153 * ElastoPlasticProblem<dim>::make_grid ()
4155 *
if (base_mesh ==
"Timoshenko beam")
4159 *
const double length = .48,
4163 * point_2(length, depth/2);
4165 * std::vector<unsigned int> repetitions(2);
4166 * repetitions[0] = 4;
4167 * repetitions[1] = 1;
4173 * give the indicators to boundaries
for specification,
4180 * |________________|
4182 * 0 to essential boundary conditions (left edge) which are as
default
4183 * 100 to the null boundaries (upper and lower edges) where we
do not need to take care of them
4184 * 5 to the natural boundaries (right edge)
for imposing the traction force
4190 *
for (; cell!=endc; ++cell)
4194 *
if ( std::fabs(cell->face(face)->center()(0)-length) < 1e-12 )
4196 * cell->face(face)->set_manifold_id(5);
4198 *
else if ( ( std::fabs(cell->face(face)->center()(1)-(depth/2)) < 1e-12 )
4200 * ( std::fabs(cell->face(face)->center()(1)-(-depth/2)) < 1e-12 ) )
4202 * cell->face(face)->set_manifold_id(100);
4208 *
triangulation.refine_global(n_initial_global_refinements);
4211 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4215 * Example 1 from the paper: Zhong Z., .... A
new numerical method
for determining
4216 * collapse load-carrying capacity of structure made of elasto-plastic material,
4217 * J. Cent. South Univ. (2014) 21: 398-404
4223 *
const double inner_radius = .1,
4224 * outer_radius = .2;
4226 * center, inner_radius, outer_radius,
4231 * give the indicators to boundaries
for specification,
4249 * 0 - inner boundary - natural boundary condition - impose the traction force
4250 * 1 - outer boundary - free boundary - we
do not need to take care of them
4251 * 2 - left boundary - essential boundary condition - constrained to move along the x direction
4252 * 3 - bottom boundary - essential boundary condition - constrained to move along the y direction
4259 *
triangulation.set_manifold (0, inner_boundary_description);
4262 *
triangulation.set_manifold (1, outer_boundary_description);
4264 *
triangulation.refine_global(n_initial_global_refinements);
4270 *
else if (base_mesh ==
"Perforated_strip_tension")
4274 * Example 2 from the paper: Zhong Z., .... A
new numerical method
for determining
4275 * collapse load-carrying capacity of structure made of elasto-plastic material,
4276 * J. Cent. South Univ. (2014) 21: 398-404
4281 *
const int dim_2d = 2;
4283 *
const double inner_radius = 0.05,
4284 * outer_radius = 0.1,
4286 * thickness = 0.004;
4299 *
const double eps = 1
e-7 * inner_radius;
4304 * center_2d, inner_radius, outer_radius,
4309 * Modify the triangulation_1
4313 * cell = triangulation_1.begin_active(),
4314 * endc = triangulation_1.end();
4315 * std::vector<bool> treated_vertices(triangulation_1.n_vertices(),
false);
4316 *
for (; cell != endc; ++cell)
4318 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4319 *
if (cell->face(f)->at_boundary() && cell->face(f)->center()(0)>eps &&
4320 * cell->face(f)->center()(1)>eps )
4324 * distance of the face center from the center
4327 *
point(0) = cell->face(f)->center()(0) - center_2d(0);
4328 *
point(1) = cell->face(f)->center()(1) - center_2d(1);
4331 *
for (
unsigned int v=0; v < GeometryInfo<dim_2d>::vertices_per_face; ++v)
4333 *
unsigned int vv = cell->face(f)->vertex_index(v);
4334 *
if (treated_vertices[vv] ==
false)
4336 * treated_vertices[vv] =
true;
4339 * cell->face(f)->vertex(v) = center_2d+
Point<dim_2d>(outer_radius,outer_radius);
4352 * Make the triangulation_2, a rectangular above the triangulation_1
4357 * point2 (outer_radius, height);
4365 * make the triangulation_2d and
refine it
4371 * Merge the two triangulation_1 and triangulation_2
4378 * Assign boundary indicators to the boundary faces
4399 * cell = triangulation_2d.begin_active(),
4400 * endc = triangulation_2d.end();
4401 *
for (; cell != endc; ++cell)
4403 *
for (
unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
4405 *
if (cell->face(f)->at_boundary())
4407 *
if ( std::fabs(cell->face(f)->center()(1)) < eps )
4409 * cell->face(f)->set_manifold_id(1);
4411 *
else if ( std::fabs(cell->face(f)->center()(0)-outer_radius) < eps )
4413 * cell->face(f)->set_manifold_id(2);
4415 *
else if ( std::fabs(cell->face(f)->center()(1)-height) < eps )
4417 * cell->face(f)->set_manifold_id(3);
4419 *
else if ( std::fabs(cell->face(f)->center()(0)) < eps )
4421 * cell->face(f)->set_manifold_id(4);
4425 * cell->face(f)->set_all_boundary_ids(10);
4435 * triangulation_2d.set_manifold (10, inner_boundary_description);
4437 * triangulation_2d.refine_global(3);
4439 * triangulation_2d.reset_manifold (10);
4444 * Extrude the triangulation_2d and make it 3
d
4454 * Assign boundary indicators to the boundary faces
4475 *
Point<dim> center(center_2d(0), center_2d(1), 0);
4480 *
for (; cell != endc; ++cell)
4482 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4484 *
if (cell->face(f)->at_boundary())
4486 * dist_vector = cell->face(f)->center() - center;
4488 *
if ( std::fabs(dist_vector[1]) < eps )
4490 * cell->face(f)->set_manifold_id(1);
4492 *
else if ( std::fabs(dist_vector[0]-outer_radius) < eps )
4494 * cell->face(f)->set_manifold_id(2);
4496 *
else if ( std::fabs(dist_vector[1]-height) < eps )
4498 * cell->face(f)->set_manifold_id(3);
4500 *
else if ( std::fabs(dist_vector[0]) < eps )
4502 * cell->face(f)->set_manifold_id(4);
4504 *
else if ( std::fabs(dist_vector[2]) < eps )
4506 * cell->face(f)->set_manifold_id(5);
4508 *
else if ( std::fabs(dist_vector[2]-thickness) < eps )
4510 * cell->face(f)->set_manifold_id(6);
4514 * cell->face(f)->set_all_boundary_ids(10);
4524 *
triangulation.set_manifold (10, inner_boundary_description);
4526 *
triangulation.refine_global(n_initial_global_refinements);
4531 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4535 * A rectangular tube made of Aluminium
4537 * approximation of beam 17250
4538 * units are in meter
4546 *
const int dim_2d = 2;
4548 *
const double length = .7,
4551 * thickness_web = 10
e-3,
4552 * thickness_flange = 10
e-3;
4560 *
const double eps = 1
e-7 * width;
4563 * Make the triangulation_b, a rectangular at the bottom of rectangular tube
4568 * point2 (width/2, -(height/2)+thickness_flange);
4570 * std::vector<unsigned int> repetitions(dim_2d);
4571 * repetitions[0] = 8;
4572 * repetitions[1] = 1;
4579 * Make the triangulation_t, a rectangular at the top of rectangular tube
4583 *
const Point<dim_2d> point1 (-width/2, (height/2)-thickness_flange),
4584 * point2 (width/2, height/2);
4586 * std::vector<unsigned int> repetitions(dim_2d);
4587 * repetitions[0] = 8;
4588 * repetitions[1] = 1;
4595 * Make the triangulation_l, a rectangular at the left of rectangular tube
4599 *
const Point<dim_2d> point1 (-width/2, -(height/2)+thickness_flange),
4600 * point2 (-(width/2)+thickness_web, (height/2)-thickness_flange);
4602 * std::vector<unsigned int> repetitions(dim_2d);
4603 * repetitions[0] = 1;
4604 * repetitions[1] = 18;
4611 * Make the triangulation_r, a rectangular at the right of rectangular tube
4615 *
const Point<dim_2d> point1 ((width/2)-thickness_web, -(height/2)+thickness_flange),
4616 * point2 (width/2, (height/2)-thickness_flange);
4618 * std::vector<unsigned int> repetitions(dim_2d);
4619 * repetitions[0] = 1;
4620 * repetitions[1] = 18;
4627 * make the triangulation_2d
4633 * merging every two triangles to make triangulation_2d
4637 * triangulation_blr;
4646 * Extrude the triangulation_2d and make it 3
d
4649 *
const unsigned int n_slices =
static_cast<int>(length*1000/20) + 1;
4655 * Assign boundary indicators to the boundary faces
4699 *
for (; cell != endc; ++cell)
4701 *
for (
unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
4703 *
if (cell->face(f)->at_boundary())
4705 * dist_vector = cell->face(f)->center() - center;
4707 *
if ( std::fabs(dist_vector[2]) < eps )
4709 * cell->face(f)->set_manifold_id(1);
4711 *
else if ( std::fabs(dist_vector[1]-(height/2)) < eps )
4713 * cell->face(f)->set_manifold_id(2);
4717 * cell->face(f)->set_all_boundary_ids(0);
4726 *
triangulation.refine_global(n_initial_global_refinements);
4734 * pcout <<
" Number of active cells: "
4744 * <a name=
"elastoplastic.cc-PlasticityContactProblemsetup_system"></a>
4745 * <h4>PlasticityContactProblem::setup_system</h4>
4749 * The next piece in the puzzle is to set up the
DoFHandler, resize
4750 * vectors and take care of various other status variables such as
4751 *
index sets and constraint matrices.
4755 * In the following, each
group of operations is put into a brace-enclosed
4756 * block that is being timed by the variable declared at the top of the
4758 * timed section, the destructor that is called at the end of the block
4762 *
template <
int dim>
4764 * ElastoPlasticProblem<dim>::setup_system ()
4770 * dof_handler.distribute_dofs(fe);
4771 * pcout <<
" Number of degrees of freedom: "
4772 * << dof_handler.n_dofs()
4775 * locally_owned_dofs = dof_handler.locally_owned_dofs();
4776 * locally_relevant_dofs =
4783 * constraints_hanging_nodes.reinit(locally_relevant_dofs);
4785 * constraints_hanging_nodes);
4786 * constraints_hanging_nodes.close();
4788 * pcout <<
" Number of active cells: "
4790 * <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
4793 * compute_dirichlet_constraints();
4799 *
if (timestep_no==1 || current_refinement_cycle!=0)
4801 * solution.reinit(locally_relevant_dofs, mpi_communicator);
4803 * incremental_displacement.reinit(locally_relevant_dofs, mpi_communicator);
4804 * newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
4805 * newton_rhs_residual.reinit(locally_owned_dofs, mpi_communicator);
4806 * fraction_of_plastic_q_points_per_cell.reinit(
triangulation.n_active_cells());
4811 * Finally, we set up sparsity patterns and matrices.
4812 * We temporarily (ab)use the system
matrix to also build the (diagonal)
4813 *
matrix that we use in eliminating degrees of freedom that are in contact
4814 * with the obstacle, but we then immediately set the Newton
matrix back
4821 * mpi_communicator);
4824 * constraints_dirichlet_and_hanging_nodes,
false,
4825 * this_mpi_process);
4827 * newton_matrix.reinit(sp);
4835 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_dirichlet_constraints"></a>
4836 * <h4>PlasticityContactProblem::compute_dirichlet_constraints</h4>
4840 * This function, broken out of the preceding one, computes the constraints
4841 * associated with Dirichlet-type boundary conditions and puts them into the
4842 * <code>constraints_dirichlet_and_hanging_nodes</code> variable by merging
4843 * with the constraints that come from hanging nodes.
4847 * As laid out in the introduction, we need to distinguish between two
4849 * - If the domain is a box, we set the displacement to zero at the bottom,
4850 * and allow vertical movement in z-direction along the sides. As
4851 * shown in the <code>make_grid()</code> function, the former corresponds
4852 * to boundary indicator 6, the latter to 8.
4853 * - If the domain is a half sphere, then we impose zero displacement along
4854 * the curved part of the boundary, associated with boundary indicator zero.
4857 *
template <
int dim>
4859 * ElastoPlasticProblem<dim>::compute_dirichlet_constraints ()
4861 * constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
4862 * constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
4864 * std::vector<bool> component_mask(dim);
4866 *
if (base_mesh ==
"Timoshenko beam")
4870 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4871 * constraints_dirichlet_and_hanging_nodes,
4874 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
4878 * the boundary x = 0
4881 * component_mask[0] =
true;
4882 * component_mask[1] =
false;
4885 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4886 * constraints_dirichlet_and_hanging_nodes,
4890 * the boundary y = 0
4893 * component_mask[0] =
false;
4894 * component_mask[1] =
true;
4897 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4898 * constraints_dirichlet_and_hanging_nodes,
4901 *
else if (base_mesh ==
"Perforated_strip_tension")
4905 * the boundary x = 0
4908 * component_mask[0] =
true;
4909 * component_mask[1] =
false;
4910 * component_mask[2] =
false;
4913 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4914 * constraints_dirichlet_and_hanging_nodes,
4918 * the boundary y = 0
4921 * component_mask[0] =
false;
4922 * component_mask[1] =
true;
4923 * component_mask[2] =
false;
4926 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4927 * constraints_dirichlet_and_hanging_nodes,
4931 * the boundary y = imposed incremental displacement
4934 * component_mask[0] =
false;
4935 * component_mask[1] =
true;
4936 * component_mask[2] =
false;
4939 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4940 * constraints_dirichlet_and_hanging_nodes,
4943 *
else if (base_mesh ==
"Cantiliver_beam_3d")
4947 * the boundary x = y = z = 0
4950 * component_mask[0] =
true;
4951 * component_mask[1] =
true;
4952 * component_mask[2] =
true;
4955 * EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
4956 * constraints_dirichlet_and_hanging_nodes,
4965 * constraints_dirichlet_and_hanging_nodes.close();
4972 * <a name=
"elastoplastic.cc-PlasticityContactProblemassemble_newton_system"></a>
4973 * <h4>PlasticityContactProblem::assemble_newton_system</h4>
4977 * Given the complexity of the problem, it may come as a bit of a surprise
4978 * that assembling the linear system we have to solve in each Newton iteration
4979 * is actually fairly straightforward. The following function builds the Newton
4980 * right hand side and Newton
matrix. It looks fairly innocent because the
4981 * heavy lifting happens in the call to
4982 * <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
4984 * constraints we have previously computed.
4987 *
template <
int dim>
4989 * ElastoPlasticProblem<dim>::
4996 *
if (base_mesh ==
"Timoshenko beam")
4998 * traction_surface_id = 5;
5000 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5002 * traction_surface_id = 0;
5004 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5006 * traction_surface_id = 2;
5020 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5021 *
const unsigned int n_q_points = quadrature_formula.size();
5022 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5025 *
const EquationData::BodyForce<dim> body_force;
5026 * std::vector<Vector<double> > body_force_values(n_q_points,
5029 *
const EquationData::
5030 * IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5031 * std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5037 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5041 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5044 * std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
5047 * cell = dof_handler.begin_active(),
5048 * endc = dof_handler.end();
5052 *
for (; cell != endc; ++cell)
5053 *
if (cell->is_locally_owned())
5055 * fe_values.reinit(cell);
5059 * fe_values[displacement].get_function_symmetric_gradients(delta_linearization_point,
5060 * incremental_strain_tensor);
5064 * For assembling the local right hand side contributions, we need
5065 * to access the prior linearized stress
value in
this quadrature
5066 *
point. To get it, we use the user pointer of
this cell that
5067 * points into the global array to the quadrature
point data
5068 * corresponding to the
first quadrature
point of the present cell,
5069 * and then add an offset corresponding to the
index of the
5070 * quadrature
point we presently consider:
5073 *
const PointHistory<dim> *local_quadrature_points_history
5074 * =
reinterpret_cast<PointHistory<dim>*
>(cell->user_pointer());
5075 *
Assert (local_quadrature_points_history >=
5076 * &quadrature_point_history.front(),
5077 * ExcInternalError());
5078 *
Assert (local_quadrature_points_history <
5079 * &quadrature_point_history.back(),
5080 * ExcInternalError());
5084 * In addition, we need the
values of the external body forces at
5085 * the quadrature points on
this cell:
5088 * body_force.vector_value_list(fe_values.get_quadrature_points(),
5089 * body_force_values);
5091 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5094 * tmp_strain_tensor_qpoint = local_quadrature_points_history[q_point].old_strain
5095 * + incremental_strain_tensor[q_point];
5099 * constitutive_law.get_linearized_stress_strain_tensors(tmp_strain_tensor_qpoint,
5100 * stress_strain_tensor_linearized,
5101 * stress_strain_tensor);
5104 *
for (
unsigned int i = 0; i < dim; ++i)
5106 * rhs_values_body_force[i] = body_force_values[q_point][i];
5109 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5113 * Having computed the stress-strain tensor and its linearization,
5114 * we can now put together the parts of the
matrix and right hand side.
5115 * In both, we need the linearized stress-strain tensor times the
5116 *
symmetric gradient of @f$\varphi_i@f$, i.e. the term @f$I_\Pi\varepsilon(\varphi_i)@f$,
5117 * so we introduce an abbreviation of
this term. Recall that the
5118 *
matrix corresponds to the bilinear form
5119 * @f$A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))@f$ in the
5120 * notation of the accompanying publication, whereas the right
5121 * hand side is @f$F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))@f$
5122 * where @f$u@f$ is the current linearization points (typically the last solution).
5123 * This might suggest that the right hand side will be zero
if the material
5124 * is completely elastic (where @f$I_\Pi=P_\Pi@f$) but
this ignores the fact
5125 * that the right hand side will also contain contributions from
5126 * non-homogeneous constraints due to the contact.
5130 * The code block that follows
this adds contributions that are due to
5131 * boundary forces, should there be any.
5135 * stress_phi_i = stress_strain_tensor_linearized
5136 * * fe_values[displacement].symmetric_gradient(i, q_point);
5138 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
5140 * * fe_values[displacement].symmetric_gradient(j, q_point)
5141 * * fe_values.JxW(q_point));
5145 * * incremental_strain_tensor[q_point] )
5147 * ( ( stress_strain_tensor
5148 * * fe_values[displacement].symmetric_gradient(i, q_point))
5149 * * tmp_strain_tensor_qpoint )
5151 * ( fe_values[displacement].
value(i, q_point)
5152 * * rhs_values_body_force )
5153 * ) * fe_values.JxW(q_point);
5158 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
5159 *
if (cell->face(face)->at_boundary()
5161 * cell->face(face)->boundary_id() == traction_surface_id)
5163 * fe_values_face.reinit(cell, face);
5165 * boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5166 * boundary_force_values);
5168 *
for (
unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
5171 *
for (
unsigned int i = 0; i < dim; ++i)
5173 * rhs_values[i] = boundary_force_values[q_point][i];
5175 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5176 * cell_rhs(i) += (fe_values_face[displacement].value(i, q_point)
5178 * * fe_values_face.JxW(q_point));
5182 * cell->get_dof_indices(local_dof_indices);
5183 * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_matrix, cell_rhs,
5184 * local_dof_indices,
5200 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_nonlinear_residual"></a>
5201 * <h4>PlasticityContactProblem::compute_nonlinear_residual</h4>
5205 * The following function computes the nonlinear residual of the equation
5206 * given the current solution (or any other linearization point). This
5207 * is needed in the linear search algorithm where we need to
try various
5208 * linear combinations of previous and current (trial) solution to
5209 * compute the (real, globalized) solution of the current Newton step.
5213 * That said, in a slight abuse of the name of the function, it actually
5214 * does significantly more. For example, it also computes the vector
5215 * that corresponds to the Newton residual but without eliminating
5216 * constrained degrees of freedom. We need
this vector to compute contact
5217 * forces and, ultimately, to compute the next active set. Likewise, by
5218 * keeping track of how many quadrature points we encounter on each cell
5219 * that show plastic yielding, we also compute the
5220 * <code>fraction_of_plastic_q_points_per_cell</code> vector that we
5221 * can later output to visualize the plastic zone. In both of these cases,
5222 * the results are not necessary as part of the line search, and so we may
5223 * be wasting a small amount of time computing them. At the same time,
this
5224 * information appears as a natural by-product of what we need to
do here
5225 * anyway, and we want to collect it once at the
end of each Newton
5226 * step, so we may as well
do it here.
5230 * The actual implementation of
this function should be rather obvious:
5233 *
template <
int dim>
5235 * ElastoPlasticProblem<dim>::
5239 *
if (base_mesh ==
"Timoshenko beam")
5241 * traction_surface_id = 5;
5243 *
else if (base_mesh ==
"Thick_tube_internal_pressure")
5245 * traction_surface_id = 0;
5247 *
else if (base_mesh ==
"Cantiliver_beam_3d")
5249 * traction_surface_id = 2;
5264 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5265 *
const unsigned int n_q_points = quadrature_formula.size();
5266 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
5268 *
const EquationData::BodyForce<dim> body_force;
5269 * std::vector<Vector<double> > body_force_values(n_q_points,
5272 *
const EquationData::
5273 * IncrementalBoundaryForce<dim> boundary_force(present_time, end_time);
5274 * std::vector<Vector<double> > boundary_force_values(n_face_q_points,
5279 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5283 * newton_rhs_residual = 0;
5285 * fraction_of_plastic_q_points_per_cell = 0;
5288 * cell = dof_handler.begin_active(),
5289 * endc = dof_handler.end();
5290 *
unsigned int cell_number = 0;
5291 *
for (; cell != endc; ++cell, ++cell_number)
5292 *
if (cell->is_locally_owned())
5294 * fe_values.reinit(cell);
5297 * std::vector<SymmetricTensor<2, dim> > strain_tensors(n_q_points);
5298 * fe_values[displacement].get_function_symmetric_gradients(linearization_point,
5301 * body_force.vector_value_list(fe_values.get_quadrature_points(),
5302 * body_force_values);
5304 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
5307 *
const bool q_point_is_plastic
5308 * = constitutive_law.get_stress_strain_tensor(strain_tensors[q_point],
5309 * stress_strain_tensor);
5310 *
if (q_point_is_plastic)
5311 * ++fraction_of_plastic_q_points_per_cell(cell_number);
5314 *
for (
unsigned int i = 0; i < dim; ++i)
5316 * rhs_values_body_force[i] = body_force_values[q_point][i];
5319 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5321 * cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5322 * * rhs_values_body_force
5324 * strain_tensors[q_point]
5325 * * stress_strain_tensor
5326 * * fe_values[displacement].symmetric_gradient(i, q_point)
5328 * * fe_values.JxW(q_point);
5332 * cell_rhs(i) += (fe_values[displacement].value(i, q_point)
5334 * * fe_values.JxW(q_point));
5338 *
for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
5339 *
if (cell->face(face)->at_boundary()
5340 * && cell->face(face)->boundary_id() == traction_surface_id)
5342 * fe_values_face.reinit(cell, face);
5344 * boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
5345 * boundary_force_values);
5347 *
for (
unsigned int q_point = 0; q_point < n_face_q_points;
5351 *
for (
unsigned int i = 0; i < dim; ++i)
5353 * rhs_values[i] = boundary_force_values[q_point][i];
5355 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5356 * cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values
5357 * * fe_values_face.JxW(q_point));
5361 * cell->get_dof_indices(local_dof_indices);
5362 * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_rhs,
5363 * local_dof_indices,
5364 * newton_rhs_residual);
5368 * fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
5380 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton_system"></a>
5381 * <h4>PlasticityContactProblem::solve_newton_system</h4>
5385 * The last piece before we can discuss the actual Newton iteration
5386 * on a single mesh is the solver
for the linear systems. There are
5387 * a couple of complications that slightly obscure the code, but
5388 * mostly it is just setup then solve. Among the complications are:
5392 * - For the hanging nodes we have to
apply
5394 * This is necessary
if a hanging node with solution
value @f$x_0@f$
5395 * has one neighbor with
value @f$x_1@f$ which is in contact with the
5396 * obstacle and one neighbor @f$x_2@f$ which is not in contact. Because
5397 * the update
for the former will be prescribed, the hanging node constraint
5398 * will have an inhomogeneity and will look like @f$x_0 = x_1/2 + \text{gap}/2@f$.
5399 * So the corresponding entries in the
5400 * ride-hang-side are non-zero with a
5403 * - Like in @ref step_40
"step-40", we need to shuffle between vectors that
do and
do
5404 *
do not have ghost elements when solving or
using the solution.
5408 * The rest of the function is similar to @ref step_40
"step-40" and
5409 * @ref step_41
"step-41" except that we use a BiCGStab solver
5410 * instead of CG. This is due to the fact that
for very small hardening
5411 * parameters @f$\gamma@f$, the linear system becomes almost semidefinite though
5412 * still
symmetric. BiCGStab appears to have an easier time with such linear
5416 *
template <
int dim>
5418 * ElastoPlasticProblem<dim>::solve_newton_system ()
5423 * distributed_solution = incremental_displacement;
5425 * constraints_hanging_nodes.set_zero(distributed_solution);
5426 * constraints_hanging_nodes.set_zero(newton_rhs);
5430 * ------- Solver Bicgstab --- Preconditioner AMG -------------------
5437 * std::vector<std::vector<bool> > constant_modes;
5447 * additional_data.
w_cycle =
false;
5454 * preconditioner.initialize(newton_matrix, additional_data);
5469 *
const double relative_accuracy = 1
e-2;
5470 *
const double solver_tolerance = relative_accuracy
5471 * * newton_matrix.residual(tmp, distributed_solution,
5477 * solver_tolerance);
5479 * solver.solve(newton_matrix, distributed_solution,
5480 * newton_rhs, preconditioner);
5484 * pcout <<
" Error: " << solver_control.initial_value()
5485 * <<
" -> " << solver_control.last_value() <<
" in "
5486 * << solver_control.last_step() <<
" Bicgstab iterations."
5492 * ------- Solver CG --- Preconditioner SSOR -------------------
5500 * preconditioner.initialize(newton_matrix, additional_data);
5510 *
const double relative_accuracy = 1
e-8;
5513 *
const double relative_accuracy = 1
e-2;
5514 *
const double solver_tolerance = relative_accuracy
5515 * * newton_matrix.residual(tmp, distributed_solution,
5521 * solver_tolerance);
5525 * solver_tolerance);
5527 * solver.solve(newton_matrix, distributed_solution,
5528 * newton_rhs, preconditioner);
5530 * pcout <<
" Error: " << solver_control.initial_value()
5531 * <<
" -> " << solver_control.last_value() <<
" in "
5532 * << solver_control.last_step() <<
" CG iterations."
5537 * ........................................................
5543 * constraints_dirichlet_and_hanging_nodes.distribute(distributed_solution);
5545 * incremental_displacement = distributed_solution;
5552 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton"></a>
5553 * <h4>PlasticityContactProblem::solve_newton</h4>
5557 * This is,
finally, the function that implements the damped Newton method
5558 * on the current mesh. There are two nested loops: the outer
loop for the Newton
5559 * iteration and the inner
loop for the line search which
5560 * will be used only
if necessary. To obtain a good and reasonable
5561 * starting
value we solve an elastic problem in very
first Newton step on each
5562 * mesh (or only on the
first mesh
if we transfer solutions between meshes). We
5563 *
do so by setting the yield stress to an unreasonably large
value in these
5564 * iterations and then setting it back to the correct
value in subsequent
5569 * Other than
this, the top part of
this function should be reasonably
5573 *
template <
int dim>
5575 * ElastoPlasticProblem<dim>::solve_newton ()
5584 *
double residual_norm;
5585 *
double previous_residual_norm = -std::numeric_limits<double>::max();
5588 * previous_disp_norm = 0;
5590 *
const double correct_sigma = sigma_0;
5592 *
const unsigned int max_newton_iter = 100;
5594 *
for (
unsigned int newton_step = 1; newton_step <= max_newton_iter; ++newton_step)
5596 *
if (newton_step == 1
5598 * ((transfer_solution && timestep_no == 1)
5600 * !transfer_solution))
5601 * constitutive_law.set_sigma_0(1e+10);
5603 * constitutive_law.set_sigma_0(correct_sigma);
5605 * pcout <<
" " << std::endl;
5606 * pcout <<
" Newton iteration " << newton_step << std::endl;
5608 * pcout <<
" Assembling system... " << std::endl;
5609 * newton_matrix = 0;
5611 * newton_rhs_residual = 0;
5613 * tmp_solution = solution;
5614 * tmp_solution += incremental_displacement;
5615 * assemble_newton_system(tmp_solution,
5616 * incremental_displacement);
5618 * pcout <<
" Solving system... " << std::endl;
5619 * solve_newton_system();
5623 * It gets a bit more hairy after we have computed the
5624 * trial solution @f$\tilde{\mathbf u}@f$ of the current Newton step.
5625 * We handle a highly nonlinear problem so we have to damp
5626 * Newton
's method using a line search. To understand how we do this,
5627 * recall that in our formulation, we compute a trial solution
5628 * in each Newton step and not the update between old and new solution.
5629 * Since the solution set is a convex set, we will use a line
5630 * search that tries linear combinations of the
5631 * previous and the trial solution to guarantee that the
5632 * damped solution is in our solution set again.
5633 * At most we apply 5 damping steps.
5637 * There are exceptions to when we use a line search. First,
5638 * if this is the first Newton step on any mesh, then we don't have
5639 * any
point to compare the residual to, so we
always accept a full
5640 * step. Likewise,
if this is the
second Newton step on the
first mesh (or
5641 * the
second on any mesh
if we don
't transfer solutions from
5642 * mesh to mesh), then we have computed the first of these steps using
5643 * just an elastic model (see how we set the yield stress sigma to
5644 * an unreasonably large value above). In this case, the first Newton
5645 * solution was a purely elastic one, the second one a plastic one,
5646 * and any linear combination would not necessarily be expected to
5647 * lie in the feasible set -- so we just accept the solution we just
5652 * In either of these two cases, we bypass the line search and just
5653 * update residual and other vectors as necessary.
5656 * if ((newton_step==1)
5658 * (transfer_solution && newton_step == 2 && current_refinement_cycle == 0)
5660 * (!transfer_solution && newton_step == 2))
5662 * tmp_solution = solution;
5663 * tmp_solution += incremental_displacement;
5664 * compute_nonlinear_residual(tmp_solution);
5665 * old_solution = incremental_displacement;
5667 * residual = newton_rhs_residual;
5669 * residual.compress(VectorOperation::insert);
5671 * residual_norm = residual.l2_norm();
5673 * pcout << " Accepting Newton solution with residual: "
5674 * << residual_norm << std::endl;
5678 * for (unsigned int i = 0; i < 5; ++i)
5680 * distributed_solution = incremental_displacement;
5682 * const double alpha = std::pow(0.5, static_cast<double>(i));
5683 * tmp_vector = old_solution;
5684 * tmp_vector.sadd(1 - alpha, alpha, distributed_solution);
5686 * TimerOutput::Scope t(computing_timer, "Residual and lambda");
5688 * locally_relevant_tmp_vector = tmp_vector;
5689 * tmp_solution = solution;
5690 * tmp_solution += locally_relevant_tmp_vector;
5691 * compute_nonlinear_residual(tmp_solution);
5692 * residual = newton_rhs_residual;
5694 * residual.compress(VectorOperation::insert);
5696 * residual_norm = residual.l2_norm();
5698 * pcout << " Residual of the system: "
5699 * << residual_norm << std::endl
5700 * << " with a damping parameter alpha = " << alpha
5703 * if (residual_norm < previous_residual_norm)
5707 * incremental_displacement = tmp_vector;
5708 * old_solution = incremental_displacement;
5711 * disp_norm = incremental_displacement.l2_norm();
5716 * The final step is to check for convergence. If the residual is
5717 * less than a threshold of @f$10^{-10}@f$, then we terminate
5718 * the iteration on the current mesh:
5719 * if (residual_norm < 1e-10)
5722 * if (residual_norm < 1e-7)
5725 * pcout << " difference of two consecutive incremental displacement l2 norm : "
5726 * << std::abs(disp_norm - previous_disp_norm) << std::endl;
5727 * if ( std::abs(disp_norm - previous_disp_norm) < 1e-10 &&
5728 * (residual_norm < 1e-5 || std::abs(residual_norm - previous_residual_norm)<1e-9) )
5730 * pcout << " Convergence by difference of two consecutive solution! " << std::endl;
5735 * previous_residual_norm = residual_norm;
5736 * previous_disp_norm = disp_norm;
5743 * <a name="elastoplastic.cc-PlasticityContactProblemcompute_error"></a>
5744 * <h4>PlasticityContactProblem::compute_error</h4>
5750 * template <int dim>
5752 * ElastoPlasticProblem<dim>::compute_error ()
5754 * TrilinosWrappers::MPI::Vector tmp_solution(locally_owned_dofs, mpi_communicator);
5755 * tmp_solution = solution;
5756 * tmp_solution += incremental_displacement;
5758 * estimated_error_per_cell.reinit (triangulation.n_active_cells());
5759 * if (error_estimation_strategy == ErrorEstimationStrategy::kelly_error)
5761 * using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
5763 * KellyErrorEstimator<dim>::estimate(dof_handler,
5764 * QGauss<dim - 1>(fe.degree + 2),
5765 * std::map<types::boundary_id, const Function<dim> *>(),
5767 * estimated_error_per_cell);
5770 * else if (error_estimation_strategy == ErrorEstimationStrategy::residual_error)
5772 * compute_error_residual(tmp_solution);
5775 * else if (error_estimation_strategy == ErrorEstimationStrategy::weighted_residual_error)
5779 * make a non-parallel copy of tmp_solution
5782 * Vector<double> copy_solution(tmp_solution);
5786 * the dual function definition (it should be defined previously, e.g. input file)
5789 * if (base_mesh == "Timoshenko beam")
5791 * double length = .48,
5794 * const Point<dim> evaluation_point(length, -depth/2);
5796 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5798 * DualSolver<dim> dual_solver(triangulation, fe,
5800 * constitutive_law, dual_functional,
5801 * timestep_no, output_dir, base_mesh,
5802 * present_time, end_time);
5804 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5807 * else if (base_mesh == "Thick_tube_internal_pressure")
5809 * const unsigned int face_id = 0;
5810 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5811 * for (unsigned int i=0; i!=dim; ++i)
5813 * comp_stress[i].resize(dim);
5814 * for (unsigned int j=0; j!=dim; ++j)
5816 * comp_stress[i][j] = 1;
5820 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5822 * DualSolver<dim> dual_solver(triangulation, fe,
5824 * constitutive_law, dual_functional,
5825 * timestep_no, output_dir, base_mesh,
5826 * present_time, end_time);
5828 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5831 * else if (base_mesh == "Perforated_strip_tension")
5835 * .........................................
5836 * Mean stress_yy over the bottom boundary
5839 * const unsigned int face_id = 1;
5840 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5841 * for (unsigned int i=0; i!=dim; ++i)
5843 * comp_stress[i].resize(dim);
5844 * for (unsigned int j=0; j!=dim; ++j)
5846 * comp_stress[i][j] = 0;
5849 * comp_stress[1][1] = 1;
5851 * DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);
5855 * .........................................
5861 * DualSolver<dim> dual_solver(triangulation, fe,
5863 * constitutive_law, dual_functional,
5864 * timestep_no, output_dir, base_mesh,
5865 * present_time, end_time);
5867 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5870 * else if (base_mesh == "Cantiliver_beam_3d")
5874 * Quantity of interest:
5875 * -----------------------------------------------------------
5876 * displacement at Point A (x=0, y=height/2, z=length)
5880 * const double length = .7,
5883 * const Point<dim> evaluation_point(0, height/2, length);
5885 * DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
5890 * -----------------------------------------------------------
5891 * Mean stress at the specified domain is of interest.
5892 * The interest domains are located on the bottom and top of the flanges
5893 * close to the clamped face, z = 0
5894 * top domain: height/2 - thickness_flange <= y <= height/2
5895 * 0 <= z <= 2 * thickness_flange
5896 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
5897 * 0 <= z <= 2 * thickness_flange
5903 * std::vector<std::vector<unsigned int> > comp_stress(dim);
5904 * for (unsigned int i=0; i!=dim; ++i)
5906 * comp_stress[i].resize(dim);
5907 * for (unsigned int j=0; j!=dim; ++j)
5909 * comp_stress[i][j] = 1;
5912 * DualFunctional::MeanStressDomain<dim> dual_functional(base_mesh, comp_stress);
5916 * -----------------------------------------------------------
5922 * DualSolver<dim> dual_solver(triangulation, fe,
5924 * constitutive_law, dual_functional,
5925 * timestep_no, output_dir, base_mesh,
5926 * present_time, end_time);
5928 * dual_solver.compute_error_DWR (estimated_error_per_cell);
5933 * AssertThrow(false, ExcNotImplemented());
5940 * AssertThrow(false, ExcNotImplemented());
5944 * relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm();
5946 * pcout << "Estimated relative error = " << relative_error << std::endl;
5950 * template <int dim>
5952 * ElastoPlasticProblem<dim>::compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution)
5954 * FEValues<dim> fe_values(fe, quadrature_formula,
5956 * update_gradients |
5958 * update_quadrature_points |
5959 * update_JxW_values);
5961 * const unsigned int n_q_points = quadrature_formula.size();
5962 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
5963 * SymmetricTensor<4, dim> stress_strain_tensor_linearized;
5964 * SymmetricTensor<4, dim> stress_strain_tensor;
5965 * Tensor<5, dim> stress_strain_tensor_grad;
5966 * std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
5967 * for (unsigned int i=0; i!=n_q_points; ++i)
5969 * cell_hessians[i].resize (dim);
5971 * const EquationData::BodyForce<dim> body_force;
5973 * std::vector<Vector<double> > body_force_values (n_q_points, Vector<double>(dim));
5974 * const FEValuesExtractors::Vector displacement(0);
5977 * FEFaceValues<dim> fe_face_values_cell(fe, face_quadrature_formula,
5979 * update_quadrature_points|
5980 * update_gradients |
5981 * update_JxW_values |
5982 * update_normal_vectors),
5983 * fe_face_values_neighbor (fe, face_quadrature_formula,
5985 * update_gradients |
5986 * update_JxW_values |
5987 * update_normal_vectors);
5988 * FESubfaceValues<dim> fe_subface_values_cell (fe, face_quadrature_formula,
5989 * update_gradients);
5991 * const unsigned int n_face_q_points = face_quadrature_formula.size();
5992 * std::vector<Vector<double> > jump_residual (n_face_q_points, Vector<double>(dim));
5993 * std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
5994 * for (unsigned int i=0; i!=n_face_q_points; ++i)
5996 * cell_grads[i].resize (dim);
5998 * std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
5999 * for (unsigned int i=0; i!=n_face_q_points; ++i)
6001 * neighbor_grads[i].resize (dim);
6003 * SymmetricTensor<2, dim> q_cell_strain_tensor;
6004 * SymmetricTensor<2, dim> q_neighbor_strain_tensor;
6005 * SymmetricTensor<4, dim> cell_stress_strain_tensor;
6006 * SymmetricTensor<4, dim> neighbor_stress_strain_tensor;
6009 * typename std::map<typename DoFHandler<dim>::face_iterator, Vector<double> >
6011 * typename DoFHandler<dim>::active_cell_iterator
6012 * cell = dof_handler.begin_active(),
6013 * endc = dof_handler.end();
6014 * for (; cell!=endc; ++cell)
6015 * if (cell->is_locally_owned())
6017 * for (unsigned int face_no=0;
6018 * face_no<GeometryInfo<dim>::faces_per_cell;
6021 * face_integrals[cell->face(face_no)].reinit (dim);
6022 * face_integrals[cell->face(face_no)] = -1e20;
6026 * std::vector<Vector<float> > error_indicators_vector;
6027 * error_indicators_vector.resize( triangulation.n_active_cells(),
6028 * Vector<float>(dim) );
6032 * ----------------- estimate_some -------------------------
6035 * cell = dof_handler.begin_active();
6036 * unsigned int present_cell = 0;
6037 * for (; cell!=endc; ++cell, ++present_cell)
6038 * if (cell->is_locally_owned())
6042 * --------------- integrate_over_cell -------------------
6045 * fe_values.reinit(cell);
6046 * body_force.vector_value_list(fe_values.get_quadrature_points(),
6047 * body_force_values);
6048 * fe_values[displacement].get_function_symmetric_gradients(tmp_solution,
6050 * fe_values.get_function_hessians(tmp_solution, cell_hessians);
6052 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
6054 * constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
6055 * stress_strain_tensor_linearized,
6056 * stress_strain_tensor);
6057 * constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
6058 * cell_hessians[q_point],
6059 * stress_strain_tensor_grad);
6061 * for (unsigned int i=0; i!=dim; ++i)
6063 * error_indicators_vector[present_cell](i) +=
6064 * body_force_values[q_point](i)*fe_values.JxW(q_point);
6065 * for (unsigned int j=0; j!=dim; ++j)
6067 * for (unsigned int k=0; k!=dim; ++k)
6069 * for (unsigned int l=0; l!=dim; ++l)
6071 * error_indicators_vector[present_cell](i) +=
6072 * ( stress_strain_tensor[i][j][k][l]*
6073 * 0.5*(cell_hessians[q_point][k][l][j]
6075 * cell_hessians[q_point][l][k][j])
6076 * + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l]
6078 * fe_values.JxW(q_point);
6088 * -------------------------------------------------------
6089 * compute face_integrals
6092 * for (unsigned int face_no=0;
6093 * face_no<GeometryInfo<dim>::faces_per_cell;
6096 * if (cell->face(face_no)->at_boundary())
6098 * for (unsigned int id=0; id!=dim; ++id)
6100 * face_integrals[cell->face(face_no)](id) = 0;
6105 * if ((cell->neighbor(face_no)->has_children() == false) &&
6106 * (cell->neighbor(face_no)->level() == cell->level()) &&
6107 * (cell->neighbor(face_no)->index() < cell->index()))
6110 * if (cell->at_boundary(face_no) == false)
6111 * if (cell->neighbor(face_no)->level() < cell->level())
6115 * if (cell->face(face_no)->has_children() == false)
6119 * ------------- integrate_over_regular_face -----------
6122 * fe_face_values_cell.reinit(cell, face_no);
6123 * fe_face_values_cell.get_function_gradients (tmp_solution,
6126 * Assert (cell->neighbor(face_no).state() == IteratorState::valid,
6127 * ExcInternalError());
6128 * const unsigned int
6129 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6130 * const typename DoFHandler<dim>::active_cell_iterator
6131 * neighbor = cell->neighbor(face_no);
6133 * fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
6134 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6137 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6139 * q_cell_strain_tensor = 0.;
6140 * q_neighbor_strain_tensor = 0.;
6141 * for (unsigned int i=0; i!=dim; ++i)
6143 * for (unsigned int j=0; j!=dim; ++j)
6145 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6146 * cell_grads[q_point][j][i] );
6147 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6148 * neighbor_grads[q_point][j][i] );
6152 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6153 * cell_stress_strain_tensor);
6154 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6155 * neighbor_stress_strain_tensor);
6157 * jump_residual[q_point] = 0.;
6158 * for (unsigned int i=0; i!=dim; ++i)
6160 * for (unsigned int j=0; j!=dim; ++j)
6162 * for (unsigned int k=0; k!=dim; ++k)
6164 * for (unsigned int l=0; l!=dim; ++l)
6166 * jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
6167 * q_cell_strain_tensor[k][l]
6169 * neighbor_stress_strain_tensor[i][j][k][l]*
6170 * q_neighbor_strain_tensor[k][l] )*
6171 * fe_face_values_cell.normal_vector(q_point)[j];
6179 * Vector<double> face_integral_vector(dim);
6180 * face_integral_vector = 0;
6181 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6183 * for (unsigned int i=0; i!=dim; ++i)
6185 * face_integral_vector(i) += jump_residual[q_point](i) *
6186 * fe_face_values_cell.JxW(q_point);
6190 * Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
6191 * ExcInternalError());
6193 * for (unsigned int i=0; i!=dim; ++i)
6195 * Assert (face_integrals[cell->face(face_no)](i) == -1e20,
6196 * ExcInternalError());
6197 * face_integrals[cell->face(face_no)](i) = face_integral_vector(i);
6203 * -----------------------------------------------------
6211 * ------------- integrate_over_irregular_face ---------
6214 * const typename DoFHandler<dim>::face_iterator
6215 * face = cell->face(face_no);
6216 * const typename DoFHandler<dim>::cell_iterator
6217 * neighbor = cell->neighbor(face_no);
6218 * Assert (neighbor.state() == IteratorState::valid,
6219 * ExcInternalError());
6220 * Assert (neighbor->has_children(),
6221 * ExcInternalError());
6223 * const unsigned int
6224 * neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
6226 * for (unsigned int subface_no=0;
6227 * subface_no<face->n_children(); ++subface_no)
6229 * const typename DoFHandler<dim>::active_cell_iterator
6230 * neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
6231 * Assert (neighbor_child->face(neighbor_neighbor) ==
6232 * cell->face(face_no)->child(subface_no),
6233 * ExcInternalError());
6235 * fe_subface_values_cell.reinit (cell, face_no, subface_no);
6236 * fe_subface_values_cell.get_function_gradients (tmp_solution,
6238 * fe_face_values_neighbor.reinit (neighbor_child,
6239 * neighbor_neighbor);
6240 * fe_face_values_neighbor.get_function_gradients (tmp_solution,
6243 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6245 * q_cell_strain_tensor = 0.;
6246 * q_neighbor_strain_tensor = 0.;
6247 * for (unsigned int i=0; i!=dim; ++i)
6249 * for (unsigned int j=0; j!=dim; ++j)
6251 * q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
6252 * cell_grads[q_point][j][i] );
6253 * q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
6254 * neighbor_grads[q_point][j][i] );
6258 * constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
6259 * cell_stress_strain_tensor);
6260 * constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
6261 * neighbor_stress_strain_tensor);
6263 * jump_residual[q_point] = 0.;
6264 * for (unsigned int i=0; i!=dim; ++i)
6266 * for (unsigned int j=0; j!=dim; ++j)
6268 * for (unsigned int k=0; k!=dim; ++k)
6270 * for (unsigned int l=0; l!=dim; ++l)
6272 * jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
6273 * q_cell_strain_tensor[k][l]
6275 * neighbor_stress_strain_tensor[i][j][k][l]*
6276 * q_neighbor_strain_tensor[k][l] )*
6277 * fe_face_values_neighbor.normal_vector(q_point)[j];
6285 * Vector<double> face_integral_vector(dim);
6286 * face_integral_vector = 0;
6287 * for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
6289 * for (unsigned int i=0; i!=dim; ++i)
6291 * face_integral_vector(i) += jump_residual[q_point](i) *
6292 * fe_face_values_neighbor.JxW(q_point);
6296 * for (unsigned int i=0; i!=dim; ++i)
6298 * face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
6303 * Vector<double> sum (dim);
6305 * for (unsigned int subface_no=0;
6306 * subface_no<face->n_children(); ++subface_no)
6308 * Assert (face_integrals.find(face->child(subface_no)) !=
6309 * face_integrals.end(),
6310 * ExcInternalError());
6311 * for (unsigned int i=0; i!=dim; ++i)
6313 * Assert (face_integrals[face->child(subface_no)](i) != -1e20,
6314 * ExcInternalError());
6315 * sum(i) += face_integrals[face->child(subface_no)](i);
6318 * for (unsigned int i=0; i!=dim; ++i)
6320 * face_integrals[face](i) = sum(i);
6326 * -----------------------------------------------------
6336 * ----------------------------------------------------------
6343 * cell = dof_handler.begin_active();
6344 * for (; cell!=endc; ++cell, ++present_cell)
6345 * if (cell->is_locally_owned())
6347 * for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
6350 * Assert(face_integrals.find(cell->face(face_no)) !=
6351 * face_integrals.end(),
6352 * ExcInternalError());
6354 * for (unsigned int id=0; id!=dim; ++id)
6356 * error_indicators_vector[present_cell](id)
6357 * -= 0.5*face_integrals[cell->face(face_no)](id);
6362 * estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();
6372 * <a name="elastoplastic.cc-PlasticityContactProblemrefine_grid"></a>
6373 * <h4>PlasticityContactProblem::refine_grid</h4>
6377 * If you've made it
this far into the deal.II tutorial, the following
6378 * function refining the mesh should not pose any challenges to you
6379 * any more. It refines the mesh, either globally or
using the Kelly
6380 * error estimator, and
if so asked also transfers the solution from
6381 * the previous to the next mesh. In the latter
case, we also need
6382 * to compute the active set and other quantities again,
for which we
6383 * need the information computed by <code>compute_nonlinear_residual()</code>.
6386 * template <int dim>
6388 * ElastoPlasticProblem<dim>::refine_grid ()
6392 * ---------------------------------------------------------------
6393 * Make a field variable
for history variables to be able to
6394 * transfer the
data to the quadrature points of the
new mesh
6399 * history_dof_handler.distribute_dofs (history_fe);
6400 * std::vector< std::vector< Vector<double> > >
6401 * history_stress_field (dim, std::vector<
Vector<double> >(dim)),
6402 * local_history_stress_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6403 * local_history_stress_fe_values (dim, std::vector<
Vector<double> >(dim));
6406 * std::vector< std::vector< Vector<double> > >
6407 * history_strain_field (dim, std::vector<
Vector<double> >(dim)),
6408 * local_history_strain_values_at_qpoints (dim, std::vector<
Vector<double> >(dim)),
6409 * local_history_strain_fe_values (dim, std::vector<
Vector<double> >(dim));
6411 *
for (
unsigned int i=0; i<dim; ++i)
6412 *
for (
unsigned int j=0; j<dim; ++j)
6414 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6415 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6416 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6418 * history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6419 * local_history_strain_values_at_qpoints[i][j].reinit(quadrature_formula.size());
6420 * local_history_strain_fe_values[i][j].reinit(history_fe.dofs_per_cell);
6423 * quadrature_formula.size());
6426 * quadrature_formula, quadrature_formula,
6427 * qpoint_to_dof_matrix);
6429 * cell = dof_handler.begin_active(),
6430 * endc = dof_handler.end(),
6431 * dg_cell = history_dof_handler.begin_active();
6432 *
for (; cell!=endc; ++cell, ++dg_cell)
6433 *
if (cell->is_locally_owned())
6435 * PointHistory<dim> *local_quadrature_points_history
6436 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6437 *
Assert (local_quadrature_points_history >=
6438 * &quadrature_point_history.front(),
6439 * ExcInternalError());
6440 *
Assert (local_quadrature_points_history <
6441 * &quadrature_point_history.back(),
6442 * ExcInternalError());
6443 *
for (
unsigned int i=0; i<dim; ++i)
6444 *
for (
unsigned int j=0; j<dim; ++j)
6446 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6448 * local_history_stress_values_at_qpoints[i][j](q)
6449 * = local_quadrature_points_history[q].old_stress[i][j];
6451 * local_history_strain_values_at_qpoints[i][j](q)
6452 * = local_quadrature_points_history[q].old_strain[i][j];
6454 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
6455 * local_history_stress_values_at_qpoints[i][j]);
6456 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
6457 * history_stress_field[i][j]);
6459 * qpoint_to_dof_matrix.vmult (local_history_strain_fe_values[i][j],
6460 * local_history_strain_values_at_qpoints[i][j]);
6461 * dg_cell->set_dof_values (local_history_strain_fe_values[i][j],
6462 * history_strain_field[i][j]);
6469 * ---------------------------------------------------------------
6473 *
if (refinement_strategy == RefinementStrategy::refine_global)
6478 *
if (cell->is_locally_owned())
6479 * cell->set_refine_flag ();
6483 *
const double refine_fraction_cells = .3,
6484 * coarsen_fraction_cells = .03;
6487 *
const double refine_fraction_cells = .1,
6488 * coarsen_fraction_cells = .3;
6496 * estimated_error_per_cell,
6497 * refine_fraction_cells, coarsen_fraction_cells);
6504 * solution_transfer.prepare_for_coarsening_and_refinement(solution);
6509 *
if (transfer_solution)
6510 * incremental_displacement_transfer.prepare_for_coarsening_and_refinement(incremental_displacement);
6513 * history_stress_field_transfer1(history_dof_handler),
6514 * history_stress_field_transfer2(history_dof_handler);
6515 * history_stress_field_transfer0.prepare_for_coarsening_and_refinement(history_stress_field[0]);
6518 * history_stress_field_transfer1.prepare_for_coarsening_and_refinement(history_stress_field[1]);
6522 * history_stress_field_transfer2.prepare_for_coarsening_and_refinement(history_stress_field[2]);
6526 * history_strain_field_transfer1(history_dof_handler),
6527 * history_strain_field_transfer2(history_dof_handler);
6528 * history_strain_field_transfer0.prepare_for_coarsening_and_refinement(history_strain_field[0]);
6531 * history_strain_field_transfer1.prepare_for_coarsening_and_refinement(history_strain_field[1]);
6535 * history_strain_field_transfer2.prepare_for_coarsening_and_refinement(history_strain_field[2]);
6539 * pcout <<
" Number of active cells: "
6544 * setup_quadrature_point_history ();
6550 * distributed_solution = solution;
6553 * solution_transfer.interpolate(distributed_solution);
6554 * solution = distributed_solution;
6556 *
if (transfer_solution)
6561 * distributed_incremental_displacement = incremental_displacement;
6564 * incremental_displacement_transfer.interpolate(distributed_incremental_displacement);
6565 * incremental_displacement = distributed_incremental_displacement;
6568 * compute_nonlinear_residual(incremental_displacement);
6575 * ---------------------------------------------------
6578 * history_dof_handler.distribute_dofs (history_fe);
6586 *
for (
unsigned int i=0; i<dim; ++i)
6587 *
for (
unsigned int j=0; j<dim; ++j)
6589 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6592 * history_stress_field_transfer0.interpolate(history_stress_field[0]);
6595 * history_stress_field_transfer1.interpolate(history_stress_field[1]);
6599 * history_stress_field_transfer2.interpolate(history_stress_field[2]);
6607 *
for (
unsigned int i=0; i<dim; ++i)
6608 *
for (
unsigned int j=0; j<dim; ++j)
6610 * history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6613 * history_strain_field_transfer0.interpolate(history_strain_field[0]);
6616 * history_strain_field_transfer1.interpolate(history_strain_field[1]);
6620 * history_strain_field_transfer2.interpolate(history_strain_field[2]);
6628 * std::vector< std::vector< Vector<double> > >
6629 * distributed_history_stress_field (dim, std::vector<
Vector<double> >(dim));
6630 *
for (
unsigned int i=0; i<dim; ++i)
6631 *
for (
unsigned int j=0; j<dim; ++j)
6633 * distributed_history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
6636 * history_stress_field_transfer0.interpolate(history_stress_field[0], distributed_history_stress_field[0]);
6639 * history_stress_field_transfer1.interpolate(history_stress_field[1], distributed_history_stress_field[1]);
6643 * history_stress_field_transfer2.interpolate(history_stress_field[2], distributed_history_stress_field[2]);
6646 * history_stress_field = distributed_history_stress_field;
6653 * std::vector< std::vector< Vector<double> > >
6654 * distributed_history_strain_field (dim, std::vector<
Vector<double> >(dim));
6655 *
for (
unsigned int i=0; i<dim; ++i)
6656 *
for (
unsigned int j=0; j<dim; ++j)
6658 * distributed_history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
6661 * history_strain_field_transfer0.interpolate(history_strain_field[0], distributed_history_strain_field[0]);
6664 * history_strain_field_transfer1.interpolate(history_strain_field[1], distributed_history_strain_field[1]);
6668 * history_strain_field_transfer2.interpolate(history_strain_field[2], distributed_history_strain_field[2]);
6671 * history_strain_field = distributed_history_strain_field;
6676 * ---------------------------------------------------------------
6677 * Transfer the history
data to the quadrature points of the
new mesh
6678 * In a
final step, we have to get the
data back from the now
6679 * interpolated global field to the quadrature points on the
6680 *
new mesh. The following code will
do that:
6687 * history_fe.dofs_per_cell);
6690 * quadrature_formula,
6691 * dof_to_qpoint_matrix);
6692 * cell = dof_handler.begin_active();
6693 * endc = dof_handler.end();
6694 * dg_cell = history_dof_handler.begin_active();
6695 *
for (; cell != endc; ++cell, ++dg_cell)
6696 *
if (cell->is_locally_owned())
6698 * PointHistory<dim> *local_quadrature_points_history
6699 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6700 *
Assert (local_quadrature_points_history >=
6701 * &quadrature_point_history.front(),
6702 * ExcInternalError());
6703 *
Assert (local_quadrature_points_history <
6704 * &quadrature_point_history.back(),
6705 * ExcInternalError());
6706 *
for (
unsigned int i=0; i<dim; ++i)
6707 *
for (
unsigned int j=0; j<dim; ++j)
6709 * dg_cell->get_dof_values (history_stress_field[i][j],
6710 * local_history_stress_fe_values[i][j]);
6711 * dof_to_qpoint_matrix.vmult (local_history_stress_values_at_qpoints[i][j],
6712 * local_history_stress_fe_values[i][j]);
6714 * dg_cell->get_dof_values (history_strain_field[i][j],
6715 * local_history_strain_fe_values[i][j]);
6716 * dof_to_qpoint_matrix.vmult (local_history_strain_values_at_qpoints[i][j],
6717 * local_history_strain_fe_values[i][j]);
6718 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6720 * local_quadrature_points_history[q].old_stress[i][j]
6721 * = local_history_stress_values_at_qpoints[i][j](q);
6723 * local_quadrature_points_history[q].old_strain[i][j]
6724 * = local_history_strain_values_at_qpoints[i][j](q);
6735 * <a name=
"elastoplastic.cc-ElastoPlasticProblemsetup_quadrature_point_history"></a>
6736 * <h4>ElastoPlasticProblem::setup_quadrature_point_history</h4>
6740 * At the beginning of our computations, we needed to set up
initial values
6741 * of the history variables, such as the existing stresses in the material,
6742 * that we store in each quadrature
point. As mentioned above, we use the
6743 * <code>user_pointer</code>
for this that is available in each cell.
6747 * To put
this into larger perspective, we note that
if we had previously
6748 * available stresses in our model (which we assume
do not exist
for the
6749 * purpose of
this program), then we would need to
interpolate the field of
6750 * preexisting stresses to the quadrature points. Likewise,
if we were to
6751 * simulate elasto-plastic materials with hardening/softening, then we would
6752 * have to store additional history variables like the present yield stress
6753 * of the accumulated plastic strains in each quadrature
6754 * points. Pre-existing hardening or weakening would then be implemented by
6755 * interpolating these variables in the present function as well.
6758 *
template <
int dim>
6759 *
void ElastoPlasticProblem<dim>::setup_quadrature_point_history ()
6763 * What we need to
do here is to
first count how many quadrature points
6764 * are within the responsibility of
this processor. This, of course,
6765 * equals the number of cells that belong to
this processor times the
6766 * number of quadrature points our quadrature formula has on each cell.
6770 * For good measure, we also set all user pointers of all cells, whether
6771 * ours of not, to the null pointer. This way,
if we ever access the user
6772 * pointer of a cell which we should not have accessed, a segmentation
6773 * fault will let us know that
this should not have happened:
6776 *
unsigned int our_cells = 0;
6780 *
if (cell->is_locally_owned())
6787 * Next, allocate as many quadrature objects as we need. Since the
6788 * <code>resize</code> function does not actually shrink the amount of
6789 * allocated memory
if the requested
new size is smaller than the old
6790 *
size, we resort to a trick to
first free all memory, and then
6791 * reallocate it: we declare an empty vector as a temporary variable and
6792 * then
swap the contents of the old vector and
this temporary
6793 * variable. This makes sure that the
6794 * <code>quadrature_point_history</code> is now really empty, and we can
6795 * let the temporary variable that now holds the previous contents of the
6796 * vector go out of scope and be destroyed. In the next step. we can then
6797 * re-allocate as many elements as we need, with the vector
6798 *
default-initializing the <code>PointHistory</code> objects, which
6799 * includes setting the stress variables to zero.
6803 * std::vector<PointHistory<dim> > tmp;
6804 * tmp.swap (quadrature_point_history);
6806 * quadrature_point_history.resize (our_cells *
6807 * quadrature_formula.size());
6811 * Finally
loop over all cells again and set the user pointers from the
6812 * cells that belong to the present processor to
point to the
first
6813 * quadrature
point objects corresponding to
this cell in the vector of
6817 *
unsigned int history_index = 0;
6821 *
if (cell->is_locally_owned())
6823 * cell->set_user_pointer (&quadrature_point_history[history_index]);
6824 * history_index += quadrature_formula.size();
6829 * At the
end,
for good measure make sure that our count of elements was
6830 * correct and that we have both used up all objects we allocated
6831 * previously, and not
point to any objects beyond the
end of the
6832 * vector. Such defensive programming strategies are
always good checks to
6833 * avoid accidental errors and to guard against future changes to
this
6834 * function that forget to update all uses of a variable at the same
6835 * time. Recall that constructs
using the <code>
Assert</code> macro are
6836 * optimized away in optimized mode, so
do not affect the
run time of
6840 *
Assert (history_index == quadrature_point_history.size(),
6841 * ExcInternalError());
6847 * <a name=
"elastoplastic.cc-ElastoPlasticProblemupdate_quadrature_point_history"></a>
6848 * <h4>ElastoPlasticProblem::update_quadrature_point_history</h4>
6852 * At the
end of each time step, we should have computed an incremental
6853 * displacement update so that the material in its
new configuration
6854 * accommodates
for the difference between the external body and boundary
6855 * forces applied during
this time step minus the forces exerted through
6856 * preexisting
internal stresses. In order to have the preexisting
6857 * stresses available at the next time step, we therefore have to update the
6858 * preexisting stresses with the stresses due to the incremental
6859 * displacement computed during the present time step. Ideally, the
6860 * resulting
sum of
internal stresses would exactly counter all external
6861 * forces. Indeed, a simple experiment can make sure that
this is so:
if we
6862 * choose boundary conditions and body forces to be time
independent, then
6863 * the forcing terms (the sum of external forces and
internal stresses)
6864 * should be exactly zero. If you make
this experiment, you will realize
6865 * from the output of the
norm of the right hand side in each time step that
6866 *
this is almost the
case: it is not exactly zero, since in the
first time
6867 * step the incremental displacement and stress updates were computed
6868 * relative to the undeformed mesh, which was then deformed. In the
second
6869 * time step, we again compute displacement and stress updates, but
this
6870 * time in the deformed mesh -- there, the resulting updates are very small
6871 * but not quite zero. This can be iterated, and in each such iteration the
6872 * residual, i.e. the
norm of the right hand side vector, is reduced;
if one
6873 * makes
this little experiment, one realizes that the
norm of
this residual
6874 * decays exponentially with the number of iterations, and after an
initial
6875 * very rapid decline is reduced by roughly a factor of about 3.5 in each
6876 * iteration (
for one testcase I looked at, other testcases, and other
6877 *
numbers of unknowns change the factor, but not the exponential decay).
6881 * In a sense,
this can then be considered as a quasi-timestepping scheme to
6882 * resolve the nonlinear problem of solving large-deformation elasticity on
6883 * a mesh that is moved along in a Lagrangian manner.
6887 * Another complication is that the existing (old) stresses are defined on
6888 * the old mesh, which we will move around after updating the stresses. If
6889 *
this mesh update involves rotations of the cell, then we need to also
6890 *
rotate the updated stress, since it was computed relative to the
6891 * coordinate system of the old cell.
6895 * Thus, what we need is the following: on each cell which the present
6896 * processor owns, we need to
extract the old stress from the
data stored
6897 * with each quadrature
point, compute the stress update, add the two
6898 * together, and then
rotate the result together with the incremental
6899 * rotation computed from the incremental displacement at the present
6900 * quadrature
point. We will detail these steps below:
6903 *
template <
int dim>
6904 *
void ElastoPlasticProblem<dim>::
6905 * update_quadrature_point_history ()
6909 * First, set up an <code>
FEValues</code>
object by which we will evaluate
6910 * the displacements and the
gradients thereof at the
6911 * quadrature points, together with a vector that will hold
this
6919 *
const unsigned int n_q_points = quadrature_formula.size();
6921 * std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
6927 * Then
loop over all cells and
do the job in the cells that belong to our
6935 * cell = dof_handler.begin_active(),
6936 * endc = dof_handler.end();
6940 *
for (; cell != endc; ++cell)
6941 *
if (cell->is_locally_owned())
6945 * Next, get a pointer to the quadrature
point history
data local to
6946 * the present cell, and, as a defensive measure, make sure that
6947 *
this pointer is within the bounds of the global array:
6950 * PointHistory<dim> *local_quadrature_points_history
6951 * =
reinterpret_cast<PointHistory<dim> *
>(cell->user_pointer());
6952 *
Assert (local_quadrature_points_history >=
6953 * &quadrature_point_history.front(),
6954 * ExcInternalError());
6955 *
Assert (local_quadrature_points_history <
6956 * &quadrature_point_history.back(),
6957 * ExcInternalError());
6961 * Then initialize the <code>
FEValues</code>
object on the present
6962 * cell, and
extract the strains of the displacement at the
6966 * fe_values.reinit (cell);
6967 * fe_values[displacement].get_function_symmetric_gradients(incremental_displacement,
6968 * incremental_strain_tensor);
6972 * Then
loop over the quadrature points of
this cell:
6975 *
for (
unsigned int q=0; q<quadrature_formula.size(); ++q)
6977 * local_quadrature_points_history[q].old_strain +=
6978 * incremental_strain_tensor[q];
6980 * constitutive_law.get_stress_strain_tensor(local_quadrature_points_history[q].old_strain,
6981 * stress_strain_tensor);
6985 * The result of these operations is then written back into
6986 * the original place:
6989 * local_quadrature_points_history[q].old_stress
6990 * = stress_strain_tensor * local_quadrature_points_history[q].old_strain;
6992 * local_quadrature_points_history[q].point
6993 * = fe_values.get_quadrature_points ()[q];
7002 * <a name=
"elastoplastic.cc-PlasticityContactProblemmove_mesh"></a>
7003 * <h4>PlasticityContactProblem::move_mesh</h4>
7007 * The remaining three
functions before we get to <code>
run()</code>
7008 * have to
do with generating output. The following one is an attempt
7009 * at showing the deformed body in its deformed configuration. To
this
7010 *
end,
this function takes a displacement vector field and moves every
7011 * vertex of the (local part) of the mesh by the previously computed
7012 * displacement. We will call
this function with the current
7013 * displacement field before we generate graphical output, and we will
7014 * call it again after generating graphical output with the
negative
7015 * displacement field to undo the changes to the mesh so made.
7019 * The function itself is pretty straightforward. All we have to
do
7020 * is keep track which vertices we have already touched, as we
7021 * encounter the same vertices multiple times as we
loop over cells.
7024 *
template <
int dim>
7026 * ElastoPlasticProblem<dim>::
7029 * std::vector<bool> vertex_touched(
triangulation.n_vertices(),
false);
7032 * dof_handler.begin_active();
7033 * cell != dof_handler.end(); ++cell)
7034 *
if (cell->is_locally_owned())
7035 *
for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7036 *
if (vertex_touched[cell->vertex_index(v)] ==
false)
7038 * vertex_touched[cell->vertex_index(v)] =
true;
7041 *
for (
unsigned int d = 0;
d < dim; ++
d)
7042 * vertex_displacement[d] = displacement(cell->vertex_dof_index(v, d));
7044 * cell->vertex(v) += vertex_displacement;
7053 * <a name=
"elastoplastic.cc-PlasticityContactProblemoutput_results"></a>
7054 * <h4>PlasticityContactProblem::output_results</h4>
7058 * Next is the function we use to actually generate graphical output. The
7059 * function is a bit tedious, but not actually particularly complicated.
7060 * It moves the mesh at the top (and moves it back at the end), then
7061 * computes the contact forces along the contact surface. We can
do
7062 * so (as shown in the accompanying paper) by taking the untreated
7063 * residual vector and identifying which degrees of freedom
7064 * correspond to those with contact by asking whether they have an
7065 * inhomogeneous constraints associated with them. As
always, we need
7066 * to be mindful that we can only write into completely distributed
7067 * vectors (i.e., vectors without ghost elements) but that when we
7068 * want to generate output, we need vectors that
do indeed have
7069 * ghost entries
for all locally relevant degrees of freedom.
7072 *
template <
int dim>
7074 * ElastoPlasticProblem<dim>::output_results (
const std::string &filename_base)
7078 * pcout <<
" Writing graphical output... " << std::flush;
7082 *
const double magnified_factor = 3;
7083 * magnified_solution *= magnified_factor;
7085 * move_mesh(magnified_solution);
7092 *
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
7094 * data_out.add_data_vector(solution,
7095 * std::vector<std::string> (dim,
"displacement"),
7099 * std::vector<std::string> solution_names;
7104 * solution_names.push_back (
"displacement");
7107 * solution_names.push_back (
"x_displacement");
7108 * solution_names.push_back (
"y_displacement");
7111 * solution_names.push_back (
"x_displacement");
7112 * solution_names.push_back (
"y_displacement");
7113 * solution_names.push_back (
"z_displacement");
7119 * data_out.add_data_vector (solution, solution_names);
7124 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
7126 * data_out.add_data_vector(subdomain,
"subdomain");
7129 * data_out.add_data_vector(fraction_of_plastic_q_points_per_cell,
7130 *
"fraction_of_plastic_q_points");
7133 * data_out.build_patches();
7137 * In the remainder of the function, we generate one VTU file on
7138 * every processor, indexed by the subdomain
id of
this processor.
7139 * On the
first processor, we then also create a <code>.pvtu</code>
7140 * file that indexes <i>all</i> of the VTU files so that the entire
7141 * set of output files can be read at once. These <code>.pvtu</code>
7142 * are used by Paraview to describe an entire
parallel computation
's
7143 * output files. We then do the same again for the competitor of
7144 * Paraview, the Visit visualization program, by creating a matching
7145 * <code>.visit</code> file.
7148 * const std::string filename =
7149 * (output_dir + filename_base + "-"
7150 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7152 * std::ofstream output_vtu((filename + ".vtu").c_str());
7153 * data_out.write_vtu(output_vtu);
7154 * pcout << output_dir + filename_base << ".pvtu" << std::endl;
7157 * if (this_mpi_process == 0)
7159 * std::vector<std::string> filenames;
7160 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7161 * filenames.push_back(filename_base + "-" +
7162 * Utilities::int_to_string(i, 4) +
7165 * std::ofstream pvtu_master_output((output_dir + filename_base + ".pvtu").c_str());
7166 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7168 * std::ofstream visit_master_output((output_dir + filename_base + ".visit").c_str());
7169 * data_out.write_pvtu_record(visit_master_output, filenames);
7173 * produce eps files for mesh illustration
7176 * std::ofstream output_eps((filename + ".eps").c_str());
7178 * grid_out.write_eps(triangulation, output_eps);
7183 * Extrapolate the stresses from Gauss point to the nodes
7186 * SymmetricTensor<2, dim> stress_at_qpoint;
7188 * FE_DGQ<dim> history_fe (1);
7189 * DoFHandler<dim> history_dof_handler (triangulation);
7190 * history_dof_handler.distribute_dofs (history_fe);
7191 * std::vector< std::vector< Vector<double> > >
7192 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7193 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7194 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7195 * for (unsigned int i=0; i<dim; ++i)
7196 * for (unsigned int j=0; j<dim; ++j)
7198 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7199 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7200 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7203 * Vector<double> VM_stress_field (history_dof_handler.n_dofs()),
7204 * local_VM_stress_values_at_qpoints (quadrature_formula.size()),
7205 * local_VM_stress_fe_values (history_fe.dofs_per_cell);
7207 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7208 * quadrature_formula.size());
7209 * FETools::compute_projection_from_quadrature_points_matrix
7211 * quadrature_formula, quadrature_formula,
7212 * qpoint_to_dof_matrix);
7214 * typename DoFHandler<dim>::active_cell_iterator
7215 * cell = dof_handler.begin_active(),
7216 * endc = dof_handler.end(),
7217 * dg_cell = history_dof_handler.begin_active();
7219 * const FEValuesExtractors::Vector displacement(0);
7221 * for (; cell!=endc; ++cell, ++dg_cell)
7222 * if (cell->is_locally_owned())
7224 * PointHistory<dim> *local_quadrature_points_history
7225 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7226 * Assert (local_quadrature_points_history >=
7227 * &quadrature_point_history.front(),
7228 * ExcInternalError());
7229 * Assert (local_quadrature_points_history <
7230 * &quadrature_point_history.back(),
7231 * ExcInternalError());
7235 * Then loop over the quadrature points of this cell:
7238 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7240 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7242 * for (unsigned int i=0; i<dim; ++i)
7243 * for (unsigned int j=i; j<dim; ++j)
7245 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7248 * local_VM_stress_values_at_qpoints(q) = Evaluation::get_von_Mises_stress(stress_at_qpoint);
7253 * for (unsigned int i=0; i<dim; ++i)
7254 * for (unsigned int j=i; j<dim; ++j)
7256 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7257 * local_history_stress_values_at_qpoints[i][j]);
7258 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7259 * history_stress_field[i][j]);
7262 * qpoint_to_dof_matrix.vmult (local_VM_stress_fe_values,
7263 * local_VM_stress_values_at_qpoints);
7264 * dg_cell->set_dof_values (local_VM_stress_fe_values,
7272 * Save stresses on nodes by nodal averaging
7273 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7274 * in order to compute stresses on nodes (by applying nodal averaging)
7275 * Therefore, each vertex has one degree of freedom
7278 * FE_Q<dim> fe_1 (1);
7279 * DoFHandler<dim> dof_handler_1 (triangulation);
7280 * dof_handler_1.distribute_dofs (fe_1);
7282 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7283 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7285 * std::vector< std::vector< Vector<double> > >
7286 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7287 * for (unsigned int i=0; i<dim; ++i)
7288 * for (unsigned int j=0; j<dim; ++j)
7290 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7293 * Vector<double> VM_stress_on_vertices (dof_handler_1.n_dofs()),
7294 * counter_on_vertices (dof_handler_1.n_dofs());
7295 * VM_stress_on_vertices = 0;
7296 * counter_on_vertices = 0;
7298 * cell = dof_handler.begin_active();
7299 * dg_cell = history_dof_handler.begin_active();
7300 * typename DoFHandler<dim>::active_cell_iterator
7301 * cell_1 = dof_handler_1.begin_active();
7302 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7303 * if (cell->is_locally_owned())
7305 * dg_cell->get_dof_values (VM_stress_field,
7306 * local_VM_stress_fe_values);
7308 * for (unsigned int i=0; i<dim; ++i)
7309 * for (unsigned int j=0; j<dim; ++j)
7311 * dg_cell->get_dof_values (history_stress_field[i][j],
7312 * local_history_stress_fe_values[i][j]);
7315 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7317 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7322 * Point<dim> point1, point2;
7323 * point1 = cell_1->vertex(v);
7324 * point2 = dg_cell->vertex(v);
7325 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7332 * counter_on_vertices (dof_1_vertex) += 1;
7334 * VM_stress_on_vertices (dof_1_vertex) += local_VM_stress_fe_values (v);
7336 * for (unsigned int i=0; i<dim; ++i)
7337 * for (unsigned int j=0; j<dim; ++j)
7339 * history_stress_on_vertices[i][j](dof_1_vertex) +=
7340 * local_history_stress_fe_values[i][j](v);
7346 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7348 * VM_stress_on_vertices(id) /= counter_on_vertices(id);
7350 * for (unsigned int i=0; i<dim; ++i)
7351 * for (unsigned int j=0; j<dim; ++j)
7353 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7359 * Save figures of stresses
7362 * if (show_stresses)
7365 * DataOut<dim> data_out;
7366 * data_out.attach_dof_handler (history_dof_handler);
7369 * data_out.add_data_vector (history_stress_field[0][0], "stress_xx");
7370 * data_out.add_data_vector (history_stress_field[1][1], "stress_yy");
7371 * data_out.add_data_vector (history_stress_field[0][1], "stress_xy");
7372 * data_out.add_data_vector (VM_stress_field, "Von_Mises_stress");
7376 * data_out.add_data_vector (history_stress_field[0][2], "stress_xz");
7377 * data_out.add_data_vector (history_stress_field[1][2], "stress_yz");
7378 * data_out.add_data_vector (history_stress_field[2][2], "stress_zz");
7381 * data_out.build_patches ();
7383 * const std::string filename_base_stress = ("stress-" + filename_base);
7385 * const std::string filename =
7386 * (output_dir + filename_base_stress + "-"
7387 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7389 * std::ofstream output_vtu((filename + ".vtu").c_str());
7390 * data_out.write_vtu(output_vtu);
7391 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7393 * if (this_mpi_process == 0)
7395 * std::vector<std::string> filenames;
7396 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7397 * filenames.push_back(filename_base_stress + "-" +
7398 * Utilities::int_to_string(i, 4) +
7401 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7402 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7404 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7405 * data_out.write_pvtu_record(visit_master_output, filenames);
7412 * DataOut<dim> data_out;
7413 * data_out.attach_dof_handler (dof_handler_1);
7416 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_xx_averaged");
7417 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_yy_averaged");
7418 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_xy_averaged");
7419 * data_out.add_data_vector (VM_stress_on_vertices, "Von_Mises_stress_averaged");
7423 * data_out.add_data_vector (history_stress_on_vertices[0][2], "stress_xz_averaged");
7424 * data_out.add_data_vector (history_stress_on_vertices[1][2], "stress_yz_averaged");
7425 * data_out.add_data_vector (history_stress_on_vertices[2][2], "stress_zz_averaged");
7428 * data_out.build_patches ();
7430 * const std::string filename_base_stress = ("averaged-stress-" + filename_base);
7432 * const std::string filename =
7433 * (output_dir + filename_base_stress + "-"
7434 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7436 * std::ofstream output_vtu((filename + ".vtu").c_str());
7437 * data_out.write_vtu(output_vtu);
7438 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7440 * if (this_mpi_process == 0)
7442 * std::vector<std::string> filenames;
7443 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7444 * filenames.push_back(filename_base_stress + "-" +
7445 * Utilities::int_to_string(i, 4) +
7448 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7449 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7451 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7452 * data_out.write_pvtu_record(visit_master_output, filenames);
7459 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7467 * magnified_solution *= -1;
7468 * move_mesh(magnified_solution);
7475 * if (base_mesh == "Timoshenko beam")
7477 * const double length = .48,
7480 * Point<dim> intersted_point(length, -depth/2);
7481 * Point<dim> vertex_displacement;
7482 * bool vertex_found = false;
7484 * for (typename DoFHandler<dim>::active_cell_iterator cell =
7485 * dof_handler.begin_active();
7486 * cell != dof_handler.end(); ++cell)
7487 * if (cell->is_locally_owned() && !vertex_found)
7488 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7489 * if ( std::fabs(cell->vertex(v)[0] - intersted_point[0])<1e-6 &&
7490 * std::fabs(cell->vertex(v)[1] - intersted_point[1])<1e-6)
7492 * vertex_found = true;
7494 * for (unsigned int d = 0; d < dim; ++d)
7495 * vertex_displacement[d] = solution(cell->vertex_dof_index(v, d));
7500 * pcout << " Number of active cells: "
7501 * << triangulation.n_global_active_cells() << std::endl
7502 * << " Number of degrees of freedom: " << dof_handler.n_dofs()
7505 * AssertThrow(vertex_found, ExcInternalError());
7506 * std::cout << "Displacement at the point (" << intersted_point[0]
7507 * << ", " << intersted_point[1] << ") is "
7508 * << "(" << vertex_displacement[0]
7509 * << ", " << vertex_displacement[1] << ").\n";
7511 * Vector<double> vertex_exact_displacement(dim);
7512 * EquationData::IncrementalBoundaryValues<dim> incremental_boundary_values(present_time, end_time);
7513 * incremental_boundary_values.vector_value (intersted_point, vertex_exact_displacement);
7515 * std::cout << "Exact displacement at the point (" << intersted_point[0]
7516 * << ", " << intersted_point[1] << ") is "
7517 * << "(" << vertex_exact_displacement[0]
7518 * << ", " << vertex_exact_displacement[1] << ").\n\n";
7521 * else if (base_mesh == "Thick_tube_internal_pressure")
7523 * const double pressure (0.6*2.4e8),
7524 * inner_radius (.1);
7527 * const double pressure (1.94e8),
7528 * inner_radius (.1);
7536 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7537 * 3d and plane strain
7540 * const double mu (e_modulus / (2 * (1 + nu)));
7542 * const Point<dim> point_A(inner_radius, 0.);
7543 * Vector<double> disp_A(dim);
7547 * make a non-parallel copy of solution
7550 * Vector<double> copy_solution(solution);
7552 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
7554 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
7556 * table_results.add_value("time step", timestep_no);
7557 * table_results.add_value("Cells", triangulation.n_global_active_cells());
7558 * table_results.add_value("DoFs", dof_handler.n_dofs());
7559 * table_results.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7560 * table_results.add_value("4*mu*u_A/(sigma_0*a)", 4*mu*disp_A(0)/(sigma_0*inner_radius));
7564 * Compute stresses in the POLAR coordinates, 1- save it on Gauss points,
7565 * 2- extrapolate them to nodes and taking their avarages (nodal avaraging)
7568 * AssertThrow (dim == 2, ExcNotImplemented());
7572 * we define a rotation matrix to be able to transform the stress
7573 * from the Cartesian coordinate to the polar coordinate
7576 * Tensor<2, dim> rotation_matrix; // [cos sin; -sin cos] , sigma_r = rot * sigma * rot^T
7578 * FEValues<dim> fe_values (fe, quadrature_formula, update_quadrature_points |
7579 * update_values | update_gradients);
7581 * const unsigned int n_q_points = quadrature_formula.size();
7583 * std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
7584 * SymmetricTensor<4, dim> stress_strain_tensor;
7585 * Tensor<2, dim> stress_at_qpoint;
7587 * FE_DGQ<dim> history_fe (1);
7588 * DoFHandler<dim> history_dof_handler (triangulation);
7589 * history_dof_handler.distribute_dofs (history_fe);
7590 * std::vector< std::vector< Vector<double> > >
7591 * history_stress_field (dim, std::vector< Vector<double> >(dim)),
7592 * local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
7593 * local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
7594 * for (unsigned int i=0; i<dim; ++i)
7595 * for (unsigned int j=0; j<dim; ++j)
7597 * history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
7598 * local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
7599 * local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
7602 * FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
7603 * quadrature_formula.size());
7604 * FETools::compute_projection_from_quadrature_points_matrix
7606 * quadrature_formula, quadrature_formula,
7607 * qpoint_to_dof_matrix);
7609 * typename DoFHandler<dim>::active_cell_iterator
7610 * cell = dof_handler.begin_active(),
7611 * endc = dof_handler.end(),
7612 * dg_cell = history_dof_handler.begin_active();
7614 * const FEValuesExtractors::Vector displacement(0);
7616 * for (; cell!=endc; ++cell, ++dg_cell)
7617 * if (cell->is_locally_owned())
7619 * PointHistory<dim> *local_quadrature_points_history
7620 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
7621 * Assert (local_quadrature_points_history >=
7622 * &quadrature_point_history.front(),
7623 * ExcInternalError());
7624 * Assert (local_quadrature_points_history <
7625 * &quadrature_point_history.back(),
7626 * ExcInternalError());
7630 * Then loop over the quadrature points of this cell:
7633 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
7635 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
7639 * transform the stress from the Cartesian coordinate to the polar coordinate
7642 * const Point<dim> point = local_quadrature_points_history[q].point;
7643 * const double theta = std::atan2(point(1),point(0));
7650 * rotation_matrix[0][0] = std::cos(theta);
7651 * rotation_matrix[0][1] = std::sin(theta);
7652 * rotation_matrix[1][0] = -std::sin(theta);
7653 * rotation_matrix[1][1] = std::cos(theta);
7657 * stress in polar coordinate
7660 * stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix);
7662 * for (unsigned int i=0; i<dim; ++i)
7663 * for (unsigned int j=i; j<dim; ++j)
7665 * local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
7671 * for (unsigned int i=0; i<dim; ++i)
7672 * for (unsigned int j=i; j<dim; ++j)
7674 * qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
7675 * local_history_stress_values_at_qpoints[i][j]);
7676 * dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
7677 * history_stress_field[i][j]);
7683 * DataOut<dim> data_out;
7684 * data_out.attach_dof_handler (history_dof_handler);
7687 * data_out.add_data_vector (history_stress_field[0][0], "stress_rr");
7688 * data_out.add_data_vector (history_stress_field[1][1], "stress_tt");
7689 * data_out.add_data_vector (history_stress_field[0][1], "stress_rt");
7691 * data_out.build_patches ();
7693 * const std::string filename_base_stress = ("stress-polar-" + filename_base);
7695 * const std::string filename =
7696 * (output_dir + filename_base_stress + "-"
7697 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7699 * std::ofstream output_vtu((filename + ".vtu").c_str());
7700 * data_out.write_vtu(output_vtu);
7701 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7703 * if (this_mpi_process == 0)
7705 * std::vector<std::string> filenames;
7706 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7707 * filenames.push_back(filename_base_stress + "-" +
7708 * Utilities::int_to_string(i, 4) +
7711 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7712 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7714 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7715 * data_out.write_pvtu_record(visit_master_output, filenames);
7723 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7724 * construct a DoFHandler object based on FE_Q with 1 degree of freedom
7725 * in order to compute stresses on nodes (by applying nodal averaging)
7726 * Therefore, each vertex has one degree of freedom
7729 * FE_Q<dim> fe_1 (1);
7730 * DoFHandler<dim> dof_handler_1 (triangulation);
7731 * dof_handler_1.distribute_dofs (fe_1);
7733 * AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
7734 * ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));
7736 * std::vector< std::vector< Vector<double> > >
7737 * history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
7738 * for (unsigned int i=0; i<dim; ++i)
7739 * for (unsigned int j=0; j<dim; ++j)
7741 * history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
7744 * Vector<double> counter_on_vertices (dof_handler_1.n_dofs());
7745 * counter_on_vertices = 0;
7747 * cell = dof_handler.begin_active();
7748 * dg_cell = history_dof_handler.begin_active();
7749 * typename DoFHandler<dim>::active_cell_iterator
7750 * cell_1 = dof_handler_1.begin_active();
7751 * for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
7752 * if (cell->is_locally_owned())
7755 * for (unsigned int i=0; i<dim; ++i)
7756 * for (unsigned int j=0; j<dim; ++j)
7758 * dg_cell->get_dof_values (history_stress_field[i][j],
7759 * local_history_stress_fe_values[i][j]);
7762 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
7764 * types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);
7769 * Point<dim> point1, point2;
7770 * point1 = cell_1->vertex(v);
7771 * point2 = dg_cell->vertex(v);
7772 * AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
7779 * counter_on_vertices (dof_1_vertex) += 1;
7781 * for (unsigned int i=0; i<dim; ++i)
7782 * for (unsigned int j=0; j<dim; ++j)
7784 * history_stress_on_vertices[i][j](dof_1_vertex) +=
7785 * local_history_stress_fe_values[i][j](v);
7791 * for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
7793 * for (unsigned int i=0; i<dim; ++i)
7794 * for (unsigned int j=0; j<dim; ++j)
7796 * history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
7802 * DataOut<dim> data_out;
7803 * data_out.attach_dof_handler (dof_handler_1);
7806 * data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_rr_averaged");
7807 * data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_tt_averaged");
7808 * data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_rt_averaged");
7810 * data_out.build_patches ();
7812 * const std::string filename_base_stress = ("averaged-stress-polar-" + filename_base);
7814 * const std::string filename =
7815 * (output_dir + filename_base_stress + "-"
7816 * + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));
7818 * std::ofstream output_vtu((filename + ".vtu").c_str());
7819 * data_out.write_vtu(output_vtu);
7820 * pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;
7822 * if (this_mpi_process == 0)
7824 * std::vector<std::string> filenames;
7825 * for (unsigned int i = 0; i < n_mpi_processes; ++i)
7826 * filenames.push_back(filename_base_stress + "-" +
7827 * Utilities::int_to_string(i, 4) +
7830 * std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
7831 * data_out.write_pvtu_record(pvtu_master_output, filenames);
7833 * std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
7834 * data_out.write_pvtu_record(visit_master_output, filenames);
7841 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7847 * if ( std::abs( (present_time/end_time)*(pressure/sigma_0) - 0.6 ) <
7848 * .501*(present_timestep/end_time)*(pressure/sigma_0) )
7853 * table_results_2: presenting the stress_rr and stress_tt on the nodes of bottom edge
7856 * const unsigned int face_id = 3;
7858 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
7860 * bool evaluation_face_found = false;
7862 * typename DoFHandler<dim>::active_cell_iterator
7863 * cell = dof_handler.begin_active(),
7864 * endc = dof_handler.end(),
7865 * cell_1 = dof_handler_1.begin_active();
7866 * for (; cell!=endc; ++cell, ++cell_1)
7867 * if (cell->is_locally_owned())
7869 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7871 * if (cell->face(face)->at_boundary()
7873 * cell->face(face)->boundary_id() == face_id)
7875 * if (!evaluation_face_found)
7877 * evaluation_face_found = true;
7881 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7883 * types::global_dof_index dof_1_vertex =
7884 * cell_1->face(face)->vertex_dof_index(v, 0);
7885 * if (!vertices_found[dof_1_vertex])
7888 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
7890 * table_results_2.add_value("x coordinate", vertex_coordinate[0]);
7891 * table_results_2.add_value("stress_rr", history_stress_on_vertices[0][0](dof_1_vertex));
7892 * table_results_2.add_value("stress_tt", history_stress_on_vertices[1][1](dof_1_vertex));
7893 * table_results_2.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7895 * vertices_found[dof_1_vertex] = true;
7904 * AssertThrow(evaluation_face_found, ExcInternalError());
7908 * table_results_3: presenting the mean stress_rr of the nodes on the inner radius
7911 * const unsigned int face_id_2 = 0;
7913 * Tensor<2, dim> stress_node,
7914 * mean_stress_polar;
7915 * mean_stress_polar = 0;
7917 * std::vector<bool> vertices_found_2 (dof_handler_1.n_dofs(), false);
7918 * unsigned int no_vertices_found = 0;
7920 * evaluation_face_found = false;
7922 * cell = dof_handler.begin_active(),
7923 * endc = dof_handler.end(),
7924 * cell_1 = dof_handler_1.begin_active();
7925 * for (; cell!=endc; ++cell, ++cell_1)
7926 * if (cell->is_locally_owned())
7928 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
7930 * if (cell->face(face)->at_boundary()
7932 * cell->face(face)->boundary_id() == face_id_2)
7934 * if (!evaluation_face_found)
7936 * evaluation_face_found = true;
7940 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
7942 * types::global_dof_index dof_1_vertex =
7943 * cell_1->face(face)->vertex_dof_index(v, 0);
7944 * if (!vertices_found_2[dof_1_vertex])
7946 * for (unsigned int ir=0; ir<dim; ++ir)
7947 * for (unsigned int ic=0; ic<dim; ++ic)
7948 * stress_node[ir][ic] = history_stress_on_vertices[ir][ic](dof_1_vertex);
7950 * mean_stress_polar += stress_node;
7952 * vertices_found_2[dof_1_vertex] = true;
7953 * ++no_vertices_found;
7962 * AssertThrow(evaluation_face_found, ExcInternalError());
7964 * mean_stress_polar /= no_vertices_found;
7966 * table_results_3.add_value("time step", timestep_no);
7967 * table_results_3.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
7968 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
7969 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
7970 * table_results_3.add_value("radius", inner_radius);
7971 * table_results_3.add_value("mean stress_rr", mean_stress_polar[0][0]);
7972 * table_results_3.add_value("mean stress_tt", mean_stress_polar[1][1]);
7979 * else if (base_mesh == "Perforated_strip_tension")
7981 * const double imposed_displacement (0.00055),
7982 * inner_radius (0.05);
7987 * const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
7988 * 3d and plane strain
7992 * table_results: Demonstrates the result of displacement at the top left corner versus imposed tension
7997 * const Point<dim> point_C(0., height);
7998 * Vector<double> disp_C(dim);
8002 * make a non-parallel copy of solution
8005 * Vector<double> copy_solution(solution);
8007 * typename Evaluation::PointValuesEvaluation<dim>::
8008 * PointValuesEvaluation point_values_evaluation(point_C);
8010 * point_values_evaluation.compute (dof_handler, copy_solution, disp_C);
8012 * table_results.add_value("time step", timestep_no);
8013 * table_results.add_value("Cells", triangulation.n_global_active_cells());
8014 * table_results.add_value("DoFs", dof_handler.n_dofs());
8015 * table_results.add_value("4*mu*u_C/(sigma_0*r)", 4*mu*disp_C(1)/(sigma_0*inner_radius));
8021 * compute average sigma_yy on the bottom edge
8024 * double stress_yy_av;
8027 * const unsigned int face_id = 1;
8029 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
8030 * unsigned int no_vertices_in_face = 0;
8032 * bool evaluation_face_found = false;
8034 * typename DoFHandler<dim>::active_cell_iterator
8035 * cell = dof_handler.begin_active(),
8036 * endc = dof_handler.end(),
8037 * cell_1 = dof_handler_1.begin_active();
8038 * for (; cell!=endc; ++cell, ++cell_1)
8039 * if (cell->is_locally_owned())
8041 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
8043 * if (cell->face(face)->at_boundary()
8045 * cell->face(face)->boundary_id() == face_id)
8047 * if (!evaluation_face_found)
8049 * evaluation_face_found = true;
8053 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
8055 * types::global_dof_index dof_1_vertex =
8056 * cell_1->face(face)->vertex_dof_index(v, 0);
8057 * if (!vertices_found[dof_1_vertex])
8059 * stress_yy_av += history_stress_on_vertices[1][1](dof_1_vertex);
8060 * ++no_vertices_in_face;
8062 * vertices_found[dof_1_vertex] = true;
8071 * AssertThrow(evaluation_face_found, ExcInternalError());
8073 * stress_yy_av /= no_vertices_in_face;
8079 * table_results_2: Demonstrate the stress_yy on the nodes of bottom edge
8083 * if ( std::abs( (stress_yy_av/sigma_0) - .91 ) < .2 )
8086 * if ( (timestep_no) % 19 == 0 )
8093 * const unsigned int face_id = 1;
8095 * std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
8097 * bool evaluation_face_found = false;
8099 * typename DoFHandler<dim>::active_cell_iterator
8100 * cell = dof_handler.begin_active(),
8101 * endc = dof_handler.end(),
8102 * cell_1 = dof_handler_1.begin_active();
8103 * for (; cell!=endc; ++cell, ++cell_1)
8104 * if (cell->is_locally_owned())
8106 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
8108 * if (cell->face(face)->at_boundary()
8110 * cell->face(face)->boundary_id() == face_id)
8112 * if (!evaluation_face_found)
8114 * evaluation_face_found = true;
8118 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
8120 * types::global_dof_index dof_1_vertex =
8121 * cell_1->face(face)->vertex_dof_index(v, 0);
8123 * const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);
8125 * if (!vertices_found[dof_1_vertex] && std::abs(vertex_coordinate[2])<1.e-8)
8127 * table_results_2.add_value("x", vertex_coordinate[0]);
8128 * table_results_2.add_value("x/r", vertex_coordinate[0]/inner_radius);
8129 * table_results_2.add_value("stress_xx/sigma_0", history_stress_on_vertices[0][0](dof_1_vertex)/sigma_0);
8130 * table_results_2.add_value("stress_yy/sigma_0", history_stress_on_vertices[1][1](dof_1_vertex)/sigma_0);
8131 * table_results_2.add_value("stress_yy_av/sigma_0", stress_yy_av/sigma_0);
8132 * table_results_2.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8134 * vertices_found[dof_1_vertex] = true;
8143 * AssertThrow(evaluation_face_found, ExcInternalError());
8149 * table_results_3: Demonstrate the Stress_mean (average tensile stress)
8150 * on the bottom edge versus epsilon_yy on the bottom left corner
8154 * double strain_yy_A = 0.;
8158 * compute strain_yy_A
8159 * Since the point A is the node on the bottom left corner,
8160 * we need to work just with one element
8164 * const Point<dim> point_A(inner_radius, 0, 0);
8166 * Vector<double> local_strain_yy_values_at_qpoints (quadrature_formula.size()),
8167 * local_strain_yy_fe_values (history_fe.dofs_per_cell);
8169 * SymmetricTensor<2, dim> strain_at_qpoint;
8171 * typename DoFHandler<dim>::active_cell_iterator
8172 * cell = dof_handler.begin_active(),
8173 * endc = dof_handler.end(),
8174 * dg_cell = history_dof_handler.begin_active();
8176 * bool cell_found = false;
8178 * for (; cell!=endc; ++cell, ++dg_cell)
8179 * if (cell->is_locally_owned() && !cell_found)
8181 * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
8182 * if ( std::fabs(cell->vertex(v)[0] - point_A[0])<1e-6 &&
8183 * std::fabs(cell->vertex(v)[1] - point_A[1])<1e-6 &&
8184 * std::fabs(cell->vertex(v)[2] - point_A[2])<1e-6)
8186 * PointHistory<dim> *local_quadrature_points_history
8187 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8188 * Assert (local_quadrature_points_history >=
8189 * &quadrature_point_history.front(),
8190 * ExcInternalError());
8191 * Assert (local_quadrature_points_history <
8192 * &quadrature_point_history.back(),
8193 * ExcInternalError());
8197 * Then loop over the quadrature points of this cell:
8200 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8202 * strain_at_qpoint = local_quadrature_points_history[q].old_strain;
8204 * local_strain_yy_values_at_qpoints(q) = strain_at_qpoint[1][1];
8207 * qpoint_to_dof_matrix.vmult (local_strain_yy_fe_values,
8208 * local_strain_yy_values_at_qpoints);
8210 * strain_yy_A = local_strain_yy_fe_values (v);
8212 * cell_found = true;
8220 * table_results_3.add_value("time step", timestep_no);
8221 * table_results_3.add_value("Cells", triangulation.n_global_active_cells());
8222 * table_results_3.add_value("DoFs", dof_handler.n_dofs());
8223 * table_results_3.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
8224 * table_results_3.add_value("mean_tensile_stress/sigma_0", stress_yy_av/sigma_0);
8225 * table_results_3.add_value("E*strain_yy-A/sigma_0", e_modulus*strain_yy_A/sigma_0);
8230 * if (std::abs(present_time-end_time) < 1.e-7)
8232 * table_results_2.set_precision("Imposed u_y", 6);
8233 * table_results_3.set_precision("Imposed u_y", 6);
8237 * else if (base_mesh == "Cantiliver_beam_3d")
8239 * const double pressure (6e6),
8245 * table_results: Demonstrates the result of displacement at the top front point, Point A
8251 * Quantity of interest:
8252 * displacement at Point A (x=0, y=height/2, z=length)
8258 * const Point<dim> point_A(0, height/2, length);
8259 * Vector<double> disp_A(dim);
8263 * make a non-parallel copy of solution
8266 * Vector<double> copy_solution(solution);
8268 * Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);
8270 * point_values_evaluation.compute (dof_handler, copy_solution, disp_A);
8272 * table_results.add_value("time step", timestep_no);
8273 * table_results.add_value("Cells", triangulation.n_global_active_cells());
8274 * table_results.add_value("DoFs", dof_handler.n_dofs());
8275 * table_results.add_value("pressure", pressure*present_time/end_time);
8276 * table_results.add_value("u_A", disp_A(1));
8282 * demonstrate the location and maximum von-Mises stress in the
8283 * specified domain close to the clamped face, z = 0
8284 * top domain: height/2 - thickness_flange <= y <= height/2
8285 * 0 <= z <= 2 * thickness_flange
8286 * bottom domain: -height/2 <= y <= -height/2 + thickness_flange
8287 * 0 <= z <= 2 * thickness_flange
8293 * double VM_stress_max (0);
8294 * Point<dim> point_max;
8296 * SymmetricTensor<2, dim> stress_at_qpoint;
8298 * typename DoFHandler<dim>::active_cell_iterator
8299 * cell = dof_handler.begin_active(),
8300 * endc = dof_handler.end();
8302 * const FEValuesExtractors::Vector displacement(0);
8304 * for (; cell!=endc; ++cell)
8305 * if (cell->is_locally_owned())
8307 * PointHistory<dim> *local_quadrature_points_history
8308 * = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
8309 * Assert (local_quadrature_points_history >=
8310 * &quadrature_point_history.front(),
8311 * ExcInternalError());
8312 * Assert (local_quadrature_points_history <
8313 * &quadrature_point_history.back(),
8314 * ExcInternalError());
8318 * Then loop over the quadrature points of this cell:
8321 * for (unsigned int q=0; q<quadrature_formula.size(); ++q)
8323 * stress_at_qpoint = local_quadrature_points_history[q].old_stress;
8325 * const double VM_stress = Evaluation::get_von_Mises_stress(stress_at_qpoint);
8326 * if (VM_stress > VM_stress_max)
8328 * VM_stress_max = VM_stress;
8329 * point_max = local_quadrature_points_history[q].point;
8335 * table_results.add_value("maximum von_Mises stress", VM_stress_max);
8336 * table_results.add_value("x", point_max[0]);
8337 * table_results.add_value("y", point_max[1]);
8338 * table_results.add_value("z", point_max[2]);
8351 * <a name="elastoplastic.cc-PlasticityContactProblemrun"></a>
8352 * <h4>PlasticityContactProblem::run</h4>
8356 * As in all other tutorial programs, the <code>run()</code> function contains
8357 * the overall logic. There is not very much to it here: in essence, it
8358 * performs the loops over all mesh refinement cycles, and within each, hands
8359 * things over to the Newton solver in <code>solve_newton()</code> on the
8360 * current mesh and calls the function that creates graphical output for
8361 * the so-computed solution. It then outputs some statistics concerning both
8362 * run times and memory consumption that has been collected over the course of
8363 * computations on this mesh.
8366 * template <int dim>
8368 * ElastoPlasticProblem<dim>::run ()
8370 * computing_timer.reset();
8373 * present_timestep = 1;
8381 * ----------------------------------------------------------------
8382 * base_mesh == "Thick_tube_internal_pressure"
8386 * const Point<dim> center(0, 0);
8387 * const double inner_radius = .1,
8388 * outer_radius = .2;
8390 * const SphericalManifold<dim> inner_boundary_description(center, inner_radius);
8391 * triangulation.set_manifold (0, inner_boundary_description);
8393 * const SphericalManifold<dim> outer_boundary_description(center, outer_radius);
8394 * triangulation.set_manifold (1, outer_boundary_description);
8398 * ----------------------------------------------------------------
8399 * base_mesh == "Perforated_strip_tension"
8403 * const double inner_radius = 0.05;
8405 * const CylinderBoundary<dim> inner_boundary_description(inner_radius, 2);
8406 * triangulation.set_manifold (10, inner_boundary_description);
8410 * ----------------------------------------------------------------
8416 * setup_quadrature_point_history ();
8418 * while (present_time < end_time)
8420 * present_time += present_timestep;
8423 * if (present_time > end_time)
8425 * present_timestep -= (present_time - end_time);
8426 * present_time = end_time;
8428 * pcout << std::endl;
8429 * pcout << "Time step " << timestep_no << " at time " << present_time
8432 * relative_error = max_relative_error * 10;
8433 * current_refinement_cycle = 0;
8440 * ------------------------ Refinement based on the relative error -------------------------------
8446 * while (relative_error >= max_relative_error)
8451 * if ( (timestep_no > 1) && (current_refinement_cycle>0) && (relative_error >= max_relative_error) )
8453 * pcout << "The relative error, " << relative_error
8454 * << " , is still more than maximum relative error, "
8455 * << max_relative_error << ", but we move to the next increment.\n";
8456 * relative_error = .1 * max_relative_error;
8459 * if (relative_error >= max_relative_error)
8461 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8462 * ++current_refinement_cycle;
8470 * ------------------------ Refinement based on the number of refinement --------------------------
8474 * bool continue_loop = true;
8475 * while (continue_loop)
8480 * if ( (timestep_no == 1) && (current_refinement_cycle < 1) )
8482 * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
8483 * ++current_refinement_cycle;
8487 * continue_loop = false;
8495 * -------------------------------------------------------------------------------------------------
8501 * solution += incremental_displacement;
8503 * update_quadrature_point_history ();
8505 * output_results((std::string("solution-") +
8506 * Utilities::int_to_string(timestep_no, 4)).c_str());
8508 * computing_timer.print_summary();
8509 * computing_timer.reset();
8511 * Utilities::System::MemoryStats stats;
8512 * Utilities::System::get_memory_stats(stats);
8513 * pcout << "Peak virtual memory used, resident in kB: " << stats.VmSize << " "
8514 * << stats.VmRSS << std::endl;
8517 * if (std::abs(present_time-end_time) < 1.e-7)
8519 * const std::string filename = (output_dir + "Results");
8521 * std::ofstream output_txt((filename + ".txt").c_str());
8523 * pcout << std::endl;
8524 * table_results.write_text(output_txt);
8525 * pcout << std::endl;
8526 * table_results_2.write_text(output_txt);
8527 * pcout << std::endl;
8528 * table_results_3.write_text(output_txt);
8529 * pcout << std::endl;
8534 * if (base_mesh == "Thick_tube_internal_pressure")
8536 * triangulation.reset_manifold (0);
8537 * triangulation.reset_manifold (1);
8539 * else if (base_mesh == "Perforated_strip_tension")
8541 * triangulation.reset_manifold (10);
8550 * <a name="elastoplastic.cc-Thecodemaincodefunction"></a>
8551 * <h3>The <code>main</code> function</h3>
8555 * There really isn't much to the <code>main()</code> function. It looks
8559 *
int main (
int argc,
char *argv[])
8561 *
using namespace dealii;
8562 *
using namespace ElastoPlastic;
8568 *
const int dim = 3;
8569 * ElastoPlasticProblem<dim>::declare_parameters(prm);
8572 * std::cerr <<
"*** Call this program as <./elastoplastic input.prm>" << std::endl;
8576 * prm.parse_input(argv[1]);
8579 * ElastoPlasticProblem<dim> problem(prm);
8583 *
catch (std::exception &exc)
8585 * std::cerr << std::endl << std::endl
8586 * <<
"----------------------------------------------------"
8588 * std::cerr <<
"Exception on processing: " << std::endl
8589 * << exc.what() << std::endl
8590 * <<
"Aborting!" << std::endl
8591 * <<
"----------------------------------------------------"
8598 * std::cerr << std::endl << std::endl
8599 * <<
"----------------------------------------------------"
8601 * std::cerr <<
"Unknown exception!" << std::endl
8602 * <<
"Aborting!" << std::endl
8603 * <<
"----------------------------------------------------"
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
#define DEAL_II_VERSION_GTE(major, minor, subminor)
#define Assert(cond, exc)
#define DeclException1(Exception1, type1, outsequence)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
@ component_is_part_of_vector
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
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)
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)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
int(&) functions(const void *v1, const void *v2)
const types::boundary_id invalid_boundary_id
::SolutionTransfer< dim, VectorType, spacedim > SolutionTransfer
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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 > &)
1.9.8