158 * <a name=
"elastoplastic.cc-Includefiles"></a>
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."));
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());
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);
289 *
cell = input.begin_active(); cell != input.end(); ++cell)
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();
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))
317 *
quad.boundary_id = cell->face(f)->boundary_id();
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)
334 *
quad.boundary_id =
bid + 1;
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,
368 * }
else if (dim == 3)
383 * -----------------------------------------------
384 *
"Perforated_strip_tension"
390 * -----------------------------------------------
398 * -----------------------------------------------
426 *
<<
"The evaluation point " <<
arg1
427 *
<<
" was not found among the vertices of the present grid.");
449 *
const unsigned int dofs_per_vertex = dof_handler.get_fe().dofs_per_vertex;
455 *
cell = dof_handler.begin_active(),
456 *
endc = dof_handler.end();
461 *
for (
unsigned int vertex=0;
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));
490 * <a name=
"elastoplastic.cc-ThecodePointHistorycodeclass"></a>
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
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
614 * --------------------
618 * --------------------
622 *
kappa (E / (3 * (1 - 2 *
nu))),
623 *
mu (E / (2 * (1 +
nu))),
626 * --------------------
652 * <a name=
"elastoplastic.cc-ConstitutiveLawget_stress_strain_tensor"></a>
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)
739 * -----------------------------------------------
740 *
"Perforated_strip_tension"
743 * -----------------------------------------------
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)
782 * <a name=
"elastoplastic.cc-ConstitutiveLawget_linearized_stress_strain_tensors"></a>
880 *
const double t[2][2] = {{
cos(angle),
sin(angle) },
881 *
{-
sin(angle),
cos(angle) }
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 } };
951 *
const double t = 1-c;
1233 * -----------------------------
TimoshenkoBeam ---------------------------------------
1778 *
template <
int dim>
1783 *
const double end_time);
1791 *
vector_value_list (
const std::vector<
Point<dim> > &points,
1801 *
template <
int dim>
1804 *
const double end_time)
1808 *
end_time (end_time),
1813 *
template <
int dim>
1820 *
ExcDimensionMismatch (
values.size(), dim));
1822 *
const double eps = 1.e-7 * height;
1845 *
template <
int 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)
1861 *
template <
int dim>
1869 *
template <
int dim>
1874 *
const double end_time);
1883 *
vector_value_list (
const std::vector<
Point<dim> > &points,
1892 *
template <
int dim>
1895 *
const double end_time)
1899 *
end_time (end_time)
1903 *
template <
int dim>
1910 *
ExcDimensionMismatch (
values.size(), dim));
1917 *
template <
int 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)
1934 * -------------------------------------------------------------------------------
1943 *
template <
int dim>
1957 *
template <
int dim>
1973 *
<<
"The evaluation point " <<
arg1
1974 *
<<
" was not found among the vertices of the present grid.");
1981 *
template <
int dim>
1989 *
template <
int dim>
1999 *
const unsigned int dofs_per_vertex =
dof_handler_dual.get_fe().dofs_per_vertex;
2005 *
for (
unsigned int vertex=0;
2011 *
for (
unsigned int id=0;
id!=dofs_per_vertex; ++id)
2022 *
template <
int dim>
2038 *
<<
"The evaluation point " <<
arg1
2039 *
<<
" was not found among the vertices of the present grid.");
2046 *
template <
int dim>
2054 *
template <
int dim>
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);
2082 *
for (; cell!=
endc; ++cell)
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);
2100 *
cell->get_dof_indices (local_dof_indices);
2101 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2115 *
template <
int dim>
2136 *
template <
int dim>
2145 *
ExcDimensionMismatch (
comp_mask.size(), dim) );
2149 *
template <
int dim>
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;
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)
2184 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2198 *
for (; cell!=
endc; ++cell)
2203 *
if (cell->face(face)->at_boundary()
2205 *
cell->face(face)->boundary_id() ==
face_id)
2211 *
fe_face_values.reinit (cell, face);
2217 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2219 *
const unsigned int
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)
2247 *
template <
int dim>
2252 *
const std::vector<std::vector<unsigned int> > &
comp_stress);
2264 *
const std::vector<std::vector<unsigned int> >
comp_stress;
2268 *
template <
int dim>
2271 *
const std::vector<std::vector<unsigned int> > &
comp_stress )
2277 *
ExcDimensionMismatch (
comp_stress.size(), dim) );
2281 *
template <
int dim>
2324 *
cell = dof_handler.begin_active();
2333 *
if (
cell_dual->face(face)->at_boundary()
2342 *
fe_face_values.reinit (cell, face);
2345 *
fe_face_values[
displacement].get_function_symmetric_gradients(solution,
2361 *
for (
unsigned int k=0;
k!=dim; ++
k)
2363 *
for (
unsigned int l=0;
l!=dim; ++
l)
2382 *
cell_dual->get_dof_indices (local_dof_indices);
2397 *
template <
int dim>
2402 *
const std::vector<std::vector<unsigned int> > &
comp_stress);
2414 *
const std::vector<std::vector<unsigned int> >
comp_stress;
2418 *
template <
int dim>
2421 *
const std::vector<std::vector<unsigned int> > &
comp_stress )
2427 *
ExcDimensionMismatch (
comp_stress.size(), dim) );
2431 *
template <
int dim>
2457 *
const double height = 200
e-3,
2472 *
std::vector<SymmetricTensor<2, dim> >
strain_tensor(n_q_points);
2492 *
cell = dof_handler.begin_active();
2498 *
const double y = cell->center()[1],
2499 *
z = cell->center()[2];
2519 *
fe_values.reinit(cell);
2522 *
fe_values[
displacement].get_function_symmetric_gradients(solution,
2538 *
for (
unsigned int k=0;
k!=dim; ++
k)
2540 *
for (
unsigned int l=0;
l!=dim; ++
l)
2558 *
cell_dual->get_dof_indices (local_dof_indices);
2573 *
template <
int dim>
2593 *
template <
int dim>
2605 *
template <
int dim>
2627 *
rhs.reinit (dof_handler.n_dofs());
2629 *
const QGauss<dim-1> face_quadrature(dof_handler.get_fe().tensor_degree()+1);
2640 *
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
2643 *
ExcDimensionMismatch (dofs_per_vertex, dim) );
2662 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
2679 *
cell = dof_handler.begin_active(),
2680 *
endc = dof_handler.end();
2687 *
if (cell->face(face)->at_boundary()
2689 *
cell->face(face)->boundary_id() ==
face_id)
2706 *
fe_face_values.reinit (cell, face);
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)
2752 *
cell->get_dof_indices (local_dof_indices);
2753 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
2775 *
template <
int dim>
2784 *
const unsigned int ×tep_no,
2788 *
const double &end_time);
2806 *
const unsigned int fe_degree;
2829 *
unsigned int timestep_no;
2844 *
const unsigned int ×tep_no,
2848 *
const double &end_time)
2852 *
solution(solution),
2853 *
fe_degree(fe.tensor_degree()),
2862 *
timestep_no (timestep_no),
2866 *
end_time (end_time)
2880 *
dof_handler.distribute_dofs(fe);
2883 *
std::cout <<
" Number of degrees of freedom in dual problem: "
2921 *
std::vector<bool> component_mask(dim);
2931 *
else if (
base_mesh ==
"Thick_tube_internal_pressure")
2935 *
the boundary
x = 0
2938 *
component_mask[0] =
true;
2939 *
component_mask[1] =
false;
2947 *
the boundary
y = 0
2950 *
component_mask[0] =
false;
2951 *
component_mask[1] =
true;
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;
2975 *
the boundary
y = 0
2978 *
component_mask[0] =
false;
2979 *
component_mask[1] =
true;
2980 *
component_mask[2] =
false;
2991 *
component_mask[0] =
false;
2992 *
component_mask[1] =
true;
2993 *
component_mask[2] =
false;
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;
3043 *
cell = dof_handler.begin_active();
3050 *
fe_values.reinit(cell);
3055 *
std::vector<SymmetricTensor<2, dim> >
strain_tensor(n_q_points);
3056 *
fe_values[
displacement].get_function_symmetric_gradients(solution,
3082 *
cell_dual->get_dof_indices(local_dof_indices);
3084 *
local_dof_indices,
3106 * +++
direct solver +++++++++
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
3155 * +++++++++++++++++++++++++++++++++++++++++++++++++
3169 *
std::ofstream output (
filename.c_str());
3191 *
data_out.build_patches ();
3192 *
data_out.write_vtk (output);
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)
3299 *
typename std::map<typename DoFHandler<dim>::face_iterator,
Vector<double> >
3304 *
for (; cell!=
endc; ++cell)
3305 *
if (cell->is_locally_owned())
3307 *
for (
unsigned int face_no=0;
3326 *
unsigned int present_cell = 0;
3327 *
for (; cell!=
endc; ++cell, ++present_cell)
3328 *
if (cell->is_locally_owned())
3335 *
fe_values.reinit(cell);
3336 *
body_force.vector_value_list(fe_values.get_quadrature_points(),
3354 *
for (
unsigned int i=0; i!=dim; ++i)
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)
3384 * -------------------------------------------------------
3388 *
for (
unsigned int face_no=0;
3392 *
if (cell->face(
face_no)->at_boundary())
3394 *
for (
unsigned int id=0;
id!=dim; ++id)
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)
3424 *
const unsigned int
3427 *
neighbor = cell->neighbor(
face_no);
3437 *
for (
unsigned int i=0; i!=dim; ++i)
3439 *
for (
unsigned int j=0;
j!=dim; ++
j)
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)
3482 *
for (
unsigned int i=0; i!=dim; ++i)
3493 *
for (
unsigned int i=0; i!=dim; ++i)
3503 * -----------------------------------------------------
3517 *
neighbor = cell->neighbor(
face_no);
3520 *
Assert (neighbor->has_children(),
3523 *
const unsigned int
3547 *
for (
unsigned int i=0; i!=dim; ++i)
3549 *
for (
unsigned int j=0;
j!=dim; ++
j)
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)
3577 *
fe_face_values_neighbor.normal_vector(
q_point)[
j];
3585 *
fe_face_values_neighbor.get_function_values (
dual_weights,
3592 *
for (
unsigned int i=0; i!=dim; ++i)
3596 *
fe_face_values_neighbor.JxW(
q_point);
3600 *
for (
unsigned int i=0; i!=dim; ++i)
3615 *
for (
unsigned int i=0; i!=dim; ++i)
3622 *
for (
unsigned int i=0; i!=dim; ++i)
3630 * -----------------------------------------------------
3640 * ----------------------------------------------------------
3648 *
for (; cell!=
endc; ++cell, ++present_cell)
3649 *
if (cell->is_locally_owned())
3658 *
for (
unsigned int id=0;
id!=dim; ++id)
3676 * <a name=
"elastoplastic.cc-ThecodePlasticityContactProblemcodeclasstemplate"></a>
3686 * a Newton method,
using Trilinos
and p4est
3691 *
We create a Newton method
for the
3712 *
template <
int dim>
3731 *
void compute_error ();
3733 *
void refine_grid ();
3805 *
const unsigned int fe_degree;
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
3887 *
unsigned int timestep_no;
3926 *
typename RefinementStrategy::value refinement_strategy;
3960 * <a name=
"elastoplastic.cc-ImplementationofthecodePlasticityContactProblemcodeclass"></a>
3966 * <a name=
"elastoplastic.cc-PlasticityContactProblemdeclare_parameters"></a>
3967 * <
h4>PlasticityContactProblem::declare_parameters</
h4>
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>
4053 *
template <
int dim>
4060 *
pcout(std::cout, this_mpi_process == 0),
4066 *
fe_degree (prm.get_integer(
"polynomial degree")),
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")),
4090 *
std::string
strat = prm.get(
"refinement strategy");
4091 *
if (
strat ==
"global")
4092 *
refinement_strategy = RefinementStrategy::refine_global;
4093 *
else if (
strat ==
"percentage")
4098 *
strat = prm.get(
"error estimation strategy");
4099 *
if (
strat ==
"kelly_error")
4101 *
else if (
strat ==
"residual_error")
4103 *
else if (
strat ==
"weighted_residual_error")
4113 *
pcout <<
" Using output directory '" <<
output_dir <<
"'" << std::endl;
4114 *
pcout <<
" FE degree " << fe_degree << std::endl;
4115 *
pcout <<
" transfer solution "
4124 * <a name=
"elastoplastic.cc-PlasticityContactProblemmake_grid"></a>
4148 *
return Point<3>(in(2), in(1), -in(0));
4151 *
template <
int dim>
4159 *
const double length = .48,
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);
4211 *
else if (
base_mesh ==
"Thick_tube_internal_pressure")
4223 *
const double inner_radius = .1,
4270 *
else if (
base_mesh ==
"Perforated_strip_tension")
4283 *
const double inner_radius = 0.05,
4299 *
const double eps = 1
e-7 * inner_radius;
4316 *
for (; cell !=
endc; ++cell)
4319 *
if (cell->face(f)->at_boundary() && cell->face(f)->center()(0)>eps &&
4320 *
cell->face(f)->center()(1)>eps )
4333 *
unsigned int vv = cell->face(f)->vertex_index(v);
4401 *
for (; cell !=
endc; ++cell)
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);
4480 *
for (; cell !=
endc; ++cell)
4484 *
if (cell->face(f)->at_boundary())
4490 *
cell->face(f)->set_manifold_id(1);
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);
4502 *
cell->face(f)->set_manifold_id(4);
4506 *
cell->face(f)->set_manifold_id(5);
4510 *
cell->face(f)->set_manifold_id(6);
4514 *
cell->face(f)->set_all_boundary_ids(10);
4531 *
else if (
base_mesh ==
"Cantiliver_beam_3d")
4548 *
const double length = .7,
4560 *
const double eps = 1
e-7 * width;
4584 *
point2 (width/2, height/2);
4649 *
const unsigned int n_slices =
static_cast<int>(length*1000/20) + 1;
4699 *
for (; cell !=
endc; ++cell)
4703 *
if (cell->face(f)->at_boundary())
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);
4734 *
pcout <<
" Number of active cells: "
4744 * <a name=
"elastoplastic.cc-PlasticityContactProblemsetup_system"></a>
4762 *
template <
int dim>
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();
4788 *
pcout <<
" Number of active cells: "
4790 *
<<
" Number of degrees of freedom: " << dof_handler.n_dofs()
4804 *
newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
4821 *
mpi_communicator);
4825 *
this_mpi_process);
4835 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_dirichlet_constraints"></a>
4857 *
template <
int dim>
4864 *
std::vector<bool> component_mask(dim);
4874 *
else if (
base_mesh ==
"Thick_tube_internal_pressure")
4878 *
the boundary
x = 0
4881 *
component_mask[0] =
true;
4882 *
component_mask[1] =
false;
4890 *
the boundary
y = 0
4893 *
component_mask[0] =
false;
4894 *
component_mask[1] =
true;
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;
4918 *
the boundary
y = 0
4921 *
component_mask[0] =
false;
4922 *
component_mask[1] =
true;
4923 *
component_mask[2] =
false;
4934 *
component_mask[0] =
false;
4935 *
component_mask[1] =
true;
4936 *
component_mask[2] =
false;
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;
4972 * <a name=
"elastoplastic.cc-PlasticityContactProblemassemble_newton_system"></a>
4987 *
template <
int dim>
5000 *
else if (
base_mesh ==
"Thick_tube_internal_pressure")
5004 *
else if (
base_mesh ==
"Cantiliver_beam_3d")
5020 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5037 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5041 * std::vector<SymmetricTensor<2, dim> >
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);
5085 *
the quadrature points
on this cell:
5088 *
body_force.vector_value_list(fe_values.get_quadrature_points(),
5104 *
for (
unsigned int i = 0; i < dim; ++i)
5109 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5138 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
5159 *
if (cell->face(face)->at_boundary()
5171 *
for (
unsigned int i = 0; i < dim; ++i)
5175 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5182 *
cell->get_dof_indices(local_dof_indices);
5184 *
local_dof_indices,
5200 * <a name=
"elastoplastic.cc-PlasticityContactProblemcompute_nonlinear_residual"></a>
5233 *
template <
int dim>
5243 *
else if (
base_mesh ==
"Thick_tube_internal_pressure")
5247 *
else if (
base_mesh ==
"Cantiliver_beam_3d")
5264 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
5279 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5288 *
cell = dof_handler.begin_active(),
5289 *
endc = dof_handler.end();
5292 *
if (cell->is_locally_owned())
5294 *
fe_values.reinit(cell);
5297 *
std::vector<SymmetricTensor<2, dim> >
strain_tensors(n_q_points);
5301 *
body_force.vector_value_list(fe_values.get_quadrature_points(),
5314 *
for (
unsigned int i = 0; i < dim; ++i)
5319 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5339 *
if (cell->face(face)->at_boundary()
5351 *
for (
unsigned int i = 0; i < dim; ++i)
5355 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
5361 *
cell->get_dof_indices(local_dof_indices);
5363 *
local_dof_indices,
5380 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton_system"></a>
5416 *
template <
int dim>
5437 * std::vector<std::vector<bool> > constant_modes;
5447 * additional_data.
w_cycle =
false;
5454 * preconditioner.initialize(
newton_matrix, additional_data);
5477 * solver_tolerance);
5484 *
pcout <<
" Error: " << solver_control.initial_value()
5485 * <<
" -> " << solver_control.last_value() <<
" in "
5486 * << solver_control.last_step() <<
" Bicgstab iterations."
5500 *
preconditioner.initialize(
newton_matrix, additional_data);
5521 * solver_tolerance);
5525 *
solver_tolerance);
5530 *
pcout <<
" Error: " << solver_control.initial_value()
5531 *
<<
" -> " << solver_control.last_value() <<
" in "
5532 *
<< solver_control.last_step() <<
" CG iterations."
5537 * ........................................................
5552 * <a name=
"elastoplastic.cc-PlasticityContactProblemsolve_newton"></a>
5573 *
template <
int dim>
5605 *
pcout <<
" " << std::endl;
5608 *
pcout <<
" Assembling system... " << std::endl;
5618 *
pcout <<
" Solving system... " << std::endl;
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
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>
6392 * ---------------------------------------------------------------
6400 *
std::vector< std::vector< Vector<double> > >
6406 *
std::vector< std::vector< Vector<double> > >
6411 *
for (
unsigned int i=0; i<dim; ++i)
6412 *
for (
unsigned int j=0;
j<dim; ++
j)
6429 *
cell = dof_handler.begin_active(),
6430 *
endc = dof_handler.end(),
6433 *
if (cell->is_locally_owned())
6443 *
for (
unsigned int i=0; i<dim; ++i)
6444 *
for (
unsigned int j=0;
j<dim; ++
j)
6469 * ---------------------------------------------------------------
6473 *
if (refinement_strategy == RefinementStrategy::refine_global)
6478 *
if (cell->is_locally_owned())
6479 *
cell->set_refine_flag ();
6539 *
pcout <<
" Number of active cells: "
6575 * ---------------------------------------------------
6586 *
for (
unsigned int i=0; i<dim; ++i)
6587 *
for (
unsigned int j=0;
j<dim; ++
j)
6607 *
for (
unsigned int i=0; i<dim; ++i)
6608 *
for (
unsigned int j=0;
j<dim; ++
j)
6628 *
std::vector< std::vector< Vector<double> > >
6630 *
for (
unsigned int i=0; i<dim; ++i)
6631 *
for (
unsigned int j=0;
j<dim; ++
j)
6653 *
std::vector< std::vector< Vector<double> > >
6655 *
for (
unsigned int i=0; i<dim; ++i)
6656 *
for (
unsigned int j=0;
j<dim; ++
j)
6676 * ---------------------------------------------------------------
6692 *
cell = dof_handler.begin_active();
6693 *
endc = dof_handler.end();
6696 *
if (cell->is_locally_owned())
6706 *
for (
unsigned int i=0; i<dim; ++i)
6707 *
for (
unsigned int j=0;
j<dim; ++
j)
6735 * <a name=
"elastoplastic.cc-ElastoPlasticProblemsetup_quadrature_point_history"></a>
6758 *
template <
int dim>
6766 * number
of quadrature points
our quadrature formula
has on each cell.
6780 *
if (cell->is_locally_owned())
6803 *
std::vector<PointHistory<dim> > tmp;
6821 *
if (cell->is_locally_owned())
6847 * <a name=
"elastoplastic.cc-ElastoPlasticProblemupdate_quadrature_point_history"></a>
6903 *
template <
int dim>
6935 *
cell = dof_handler.begin_active(),
6936 *
endc = dof_handler.end();
6940 *
for (; cell !=
endc; ++cell)
6941 *
if (cell->is_locally_owned())
6966 *
fe_values.reinit (cell);
6993 *
= fe_values.get_quadrature_points ()[
q];
7002 * <a name=
"elastoplastic.cc-PlasticityContactProblemmove_mesh"></a>
7024 *
template <
int dim>
7032 *
dof_handler.begin_active();
7033 *
cell != dof_handler.end(); ++cell)
7034 *
if (cell->is_locally_owned())
7041 *
for (
unsigned int d = 0;
d < dim; ++
d)
7053 * <a name=
"elastoplastic.cc-PlasticityContactProblemoutput_results"></a>
7072 *
template <
int dim>
7078 *
pcout <<
" Writing graphical output... " << std::flush;
7089 *
data_out.attach_dof_handler(dof_handler);
7092 *
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
7094 *
data_out.add_data_vector(solution,
7095 *
std::vector<std::string> (dim,
"displacement"),
7124 *
for (
unsigned int i = 0; i <
subdomain.size(); ++i)
7126 *
data_out.add_data_vector(
subdomain,
"subdomain");
7130 *
"fraction_of_plastic_q_points");
7133 *
data_out.build_patches();
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>
8561 *
using namespace dealii;
8568 *
const int dim = 3;
8572 *
std::cerr <<
"*** Call this program as <./elastoplastic input.prm>" << std::endl;
8576 *
prm.parse_input(
argv[1]);
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
unsigned int depth_console(const unsigned int n)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define DEAL_II_VERSION_GTE(major, minor, subminor)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
@ component_is_part_of_vector
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
constexpr types::blas_int zero
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
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)
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)
constexpr types::boundary_id invalid_boundary_id
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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 > &)