From: Shahram Ghorashi Date: Thu, 26 May 2016 23:19:31 +0000 (-0500) Subject: Add initial version of the project. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f1eb7b888f93032f31a2b1ea4e34e3d4091ce74f;p=code-gallery.git Add initial version of the project. These are the files as provided by Shahram Ghorashi. --- diff --git a/goal_oriented_elastoplasticity/CMakeLists.txt b/goal_oriented_elastoplasticity/CMakeLists.txt new file mode 100644 index 0000000..372c562 --- /dev/null +++ b/goal_oriented_elastoplasticity/CMakeLists.txt @@ -0,0 +1,50 @@ +## +# CMake script for the step-42 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "elastoplastic") + +# Declare all source files the target consists of: +SET(TARGET_SRC + ${TARGET}.cc + # You can specify additional files here! + ) + +# Define the output that should be cleaned: +SET(CLEAN_UP_FILES *.vtu *.pvtu *.visit) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# +# Are all dependencies fullfilled? +# +IF( NOT DEAL_II_WITH_MPI OR + NOT DEAL_II_WITH_P4EST OR + NOT DEAL_II_WITH_TRILINOS ) + MESSAGE(FATAL_ERROR " +Error! The deal.II library found at ${DEAL_II_PATH} was not configured with + DEAL_II_WITH_MPI = ON + DEAL_II_WITH_P4EST = ON + DEAL_II_WITH_TRILINOS = ON +One or all of these are OFF in your installation but are required for this tutorial step." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) + +DEAL_II_INVOKE_AUTOPILOT() diff --git a/goal_oriented_elastoplasticity/Cantiliver_II_beam_3d.prm b/goal_oriented_elastoplasticity/Cantiliver_II_beam_3d.prm new file mode 100644 index 0000000..af9d3e8 --- /dev/null +++ b/goal_oriented_elastoplasticity/Cantiliver_II_beam_3d.prm @@ -0,0 +1,21 @@ +set polynomial degree = 1 +set number of initial refinements = 0 +set refinement strategy = percentage +set error estimation strategy = weighted_residual_error +set maximum relative error = 3.9e7 +set output directory = Results/DWR_VM_1 +set transfer solution = true +set base mesh = Cantiliver_beam_3d +set elasticity modulus = 7.e10 +set Poissons ratio = 0.3 +set yield stress = 2.43e8 +set isotropic hardening parameter = 0.0401097 +set show stresses = true + + + + +# refinement strategy : global / percentage +# error estimation strategy: kelly_error / residual_error / weighted_residual_error +# base mesh : Timoshenko beam / Thick_tube_internal_pressure +# / Perforated_strip_tension / Cantiliver_beam_3d diff --git a/goal_oriented_elastoplasticity/Thick_tube_internal_pressure.prm b/goal_oriented_elastoplasticity/Thick_tube_internal_pressure.prm new file mode 100644 index 0000000..ab42ced --- /dev/null +++ b/goal_oriented_elastoplasticity/Thick_tube_internal_pressure.prm @@ -0,0 +1,20 @@ +set polynomial degree = 1 +set number of initial refinements = 2 +set refinement strategy = percentage +set error estimation strategy = weighted_residual_error +set maximum relative error = 8e4 +set output directory = p1_adaptive +set transfer solution = true +set base mesh = Thick_tube_internal_pressure +set elasticity modulus = 2.1e11 +set Poissons ratio = 0.3 +set yield stress = 2.4e8 +set isotropic hardening parameter = 0 +set show stresses = false + + + + +# refinement strategy : global / percentage +# error estimation strategy: kelly_error / residual_error / weighted_residual_error +# base mesh : Timoshenko beam / Thick_tube_internal_pressure diff --git a/goal_oriented_elastoplasticity/doc/author b/goal_oriented_elastoplasticity/doc/author new file mode 100644 index 0000000..8f2a2a6 --- /dev/null +++ b/goal_oriented_elastoplasticity/doc/author @@ -0,0 +1 @@ +Seyed Shahram Ghorashi diff --git a/goal_oriented_elastoplasticity/doc/builds-on b/goal_oriented_elastoplasticity/doc/builds-on new file mode 100644 index 0000000..eb8f57b --- /dev/null +++ b/goal_oriented_elastoplasticity/doc/builds-on @@ -0,0 +1 @@ +step-14 step-42 diff --git a/goal_oriented_elastoplasticity/doc/entry-name b/goal_oriented_elastoplasticity/doc/entry-name new file mode 100644 index 0000000..97fbece --- /dev/null +++ b/goal_oriented_elastoplasticity/doc/entry-name @@ -0,0 +1 @@ +Goal-oriented mesh adaptivity in elastoplasticity problems diff --git a/goal_oriented_elastoplasticity/doc/tooltip b/goal_oriented_elastoplasticity/doc/tooltip new file mode 100644 index 0000000..0a56b1d --- /dev/null +++ b/goal_oriented_elastoplasticity/doc/tooltip @@ -0,0 +1 @@ +Solving 2d/3d elastoplasticity problem with linear isotropic hardening and adapting the mesh based on goal-oriented error estimation diff --git a/goal_oriented_elastoplasticity/elastoplastic.cc b/goal_oriented_elastoplasticity/elastoplastic.cc new file mode 100644 index 0000000..4b0d9de --- /dev/null +++ b/goal_oriented_elastoplasticity/elastoplastic.cc @@ -0,0 +1,7232 @@ +/* --------------------------------------------------------------------- + * $Id: elastoplastic.cc 31592 2013-11-08 16:47:28Z Ghorashi $ + * + * Copyright (C) 2012 - 2013 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * + * Authors: Seyed Shahram Ghorashi, Bauhaus-Universit\"at Weimar, 2014 + * Joerg Frohne, Texas A&M University and + * University of Siegen, 2012, 2013 + * Wolfgang Bangerth, Texas A&M University, 2012, 2013 + * Timo Heister, Texas A&M University, 2013 + */ + +// @sect3{Include files} +// The set of include files is not much of a surprise any more at this time: +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +// And here the only two new things among the header files: an include file in +// which symmetric tensors of rank 2 and 4 are implemented, as introduced in +// the introduction: +#include + +// And a header that implements filters for iterators looping over all +// cells. We will use this when selecting only those cells for output that are +// owned by the present process in a %parallel program: +#include + +#include +#include + +// This final include file provides the mkdir function +// that we will use to create a directory for output files, if necessary: +#include + +namespace ElastoPlastic +{ + using namespace dealii; + + void + extrude_triangulation(const Triangulation<2, 2> &input, + const unsigned int n_slices, + const double height, + Triangulation<3,3> &result) + { + // Assert (input.n_levels() == 1, + // ExcMessage ("The input triangulations must be coarse meshes.")); + Assert(result.n_cells()==0, ExcMessage("resultin Triangulation need to be empty upon calling extrude_triangulation.")); + Assert(height>0, ExcMessage("The height in extrude_triangulation needs to be positive.")); + Assert(n_slices>=2, ExcMessage("The number of slices in extrude_triangulation needs to be at least 2.")); + + std::vector > points(n_slices*input.n_used_vertices()); + std::vector > cells; + cells.reserve((n_slices-1)*input.n_active_cells()); + + for (unsigned int slice=0; slice &v = input.get_vertices()[i]; + points[i+slice*input.n_vertices()](0) = v(0); + points[i+slice*input.n_vertices()](1) = v(1); + points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1); + } + } + } + + for (Triangulation<2,2>::cell_iterator + cell = input.begin_active(); cell != input.end(); ++cell) + { + for (unsigned int slice=0; slice this_cell; + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + this_cell.vertices[v] + = cell->vertex_index(v)+slice*input.n_used_vertices(); + this_cell.vertices[v+GeometryInfo<2>::vertices_per_cell] + = cell->vertex_index(v)+(slice+1)*input.n_used_vertices(); + } + + this_cell.material_id = cell->material_id(); + cells.push_back(this_cell); + } + } + + SubCellData s; + types::boundary_id bid=0; + s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2); + for (Triangulation<2,2>::cell_iterator + cell = input.begin_active(); cell != input.end(); ++cell) + { + CellData<2> quad; + for (unsigned int f=0; f<4; ++f) + if (cell->at_boundary(f)) + { + quad.boundary_id = cell->face(f)->boundary_indicator(); + bid = std::max(bid, quad.boundary_id); + for (unsigned int slice=0; sliceface(f)->vertex_index(0)+slice*input.n_used_vertices(); + quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_used_vertices(); + quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_used_vertices(); + quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_used_vertices(); + s.boundary_quads.push_back(quad); + } + } + } + + for (Triangulation<2,2>::cell_iterator + cell = input.begin_active(); cell != input.end(); ++cell) + { + CellData<2> quad; + quad.boundary_id = bid + 1; + quad.vertices[0] = cell->vertex_index(0); + quad.vertices[1] = cell->vertex_index(1); + quad.vertices[2] = cell->vertex_index(2); + quad.vertices[3] = cell->vertex_index(3); + s.boundary_quads.push_back(quad); + + quad.boundary_id = bid + 2; + for (int i=0; i<4; ++i) + quad.vertices[i] += (n_slices-1)*input.n_used_vertices(); + s.boundary_quads.push_back(quad); + } + + result.create_triangulation (points, + cells, + s); + } + + namespace Evaluation + { + + + template + double get_von_Mises_stress(const SymmetricTensor<2, dim> &stress) + { + + // if (dim == 2) + // { + // von_Mises_stress = std::sqrt( stress[0][0]*stress[0][0] + // + stress[1][1]*stress[1][1] + // - stress[0][0]*stress[1][1] + // + 3*stress[0][1]*stress[0][1]); + // }else if (dim == 3) + // { + // von_Mises_stress = std::sqrt( stress[0][0]*stress[0][0] + // + stress[1][1]*stress[1][1] + // + stress[2][2]*stress[2][2] + // - stress[0][0]*stress[1][1] + // - stress[1][1]*stress[2][2] + // - stress[0][0]*stress[2][2] + // + 3*( stress[0][1]*stress[0][1] + // +stress[1][2]*stress[1][2] + // +stress[0][2]*stress[0][2]) ); + // } + + // ----------------------------------------------- + // "Perforated_strip_tension" + // plane stress +// const double von_Mises_stress = std::sqrt( stress[0][0]*stress[0][0] +// + stress[1][1]*stress[1][1] +// - stress[0][0]*stress[1][1] +// + 3*stress[0][1]*stress[0][1]); + // ----------------------------------------------- + // otherwise + // plane strain / 3d case + const double von_Mises_stress = std::sqrt(1.5) * (deviator(stress)).norm(); + // ----------------------------------------------- + + + + return von_Mises_stress; + } + + + template + class PointValuesEvaluation + { + public: + PointValuesEvaluation (const Point &evaluation_point); + + void compute (const DoFHandler &dof_handler, + const Vector &solution, + Vector &point_values); + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + private: + const Point evaluation_point; + }; + + + template + PointValuesEvaluation:: + PointValuesEvaluation (const Point &evaluation_point) + : + evaluation_point (evaluation_point) + {} + + + + template + void + PointValuesEvaluation:: + compute (const DoFHandler &dof_handler, + const Vector &solution, + Vector &point_values) + { + const unsigned int dofs_per_vertex = dof_handler.get_fe().dofs_per_vertex; + AssertThrow (point_values.size() == dofs_per_vertex, + ExcDimensionMismatch (point_values.size(), dofs_per_vertex)); + point_values = 1e20; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + bool evaluation_point_found = false; + for (; (cell!=endc) && !evaluation_point_found; ++cell) + { + if (cell->is_locally_owned() && !evaluation_point_found) + for (unsigned int vertex=0; + vertex::vertices_per_cell; + ++vertex) + { + if (cell->vertex(vertex).distance (evaluation_point) + < + cell->diameter() * 1e-8) + { + for (unsigned int id=0; id!=dofs_per_vertex; ++id) + { + point_values[id] = solution(cell->vertex_dof_index(vertex,id)); + } + + evaluation_point_found = true; + break; + } + } + } + + AssertThrow (evaluation_point_found, + ExcEvaluationPointNotFound(evaluation_point)); + } + + + } + + // @sect3{The PointHistory class} + + // As was mentioned in the introduction, we have to store the old stress in + // quadrature point so that we can compute the residual forces at this point + // during the next time step. This alone would not warrant a structure with + // only one member, but in more complicated applications, we would have to + // store more information in quadrature points as well, such as the history + // variables of plasticity, etc. In essence, we have to store everything + // that affects the present state of the material here, which in plasticity + // is determined by the deformation history variables. + // + // We will not give this class any meaningful functionality beyond being + // able to store data, i.e. there are no constructors, destructors, or other + // member functions. In such cases of `dumb' classes, we usually opt to + // declare them as struct rather than class, to + // indicate that they are closer to C-style structures than C++-style + // classes. + template + struct PointHistory + { + SymmetricTensor<2,dim> old_stress; + SymmetricTensor<2,dim> old_strain; + Point point; + }; + + + // @sect3{The ConstitutiveLaw class template} + + // This class provides an interface for a constitutive law, i.e., for the + // relationship between strain $\varepsilon(\mathbf u)$ and stress + // $\sigma$. In this example we are using an elastoplastic material behavior + // with linear, isotropic hardening. Such materials are characterized by + // Young's modulus $E$, Poisson's ratio $\nu$, the initial yield stress + // $\sigma_0$ and the isotropic hardening parameter $\gamma$. For $\gamma = + // 0$ we obtain perfect elastoplastic behavior. + // + // As explained in the paper that describes this program, the first Newton + // steps are solved with a completely elastic material model to avoid having + // to deal with both nonlinearities (plasticity and contact) at once. To this + // end, this class has a function set_sigma_0() that we use later + // on to simply set $\sigma_0$ to a very large value -- essentially + // guaranteeing that the actual stress will not exceed it, and thereby + // producing an elastic material. When we are ready to use a plastic model, we + // set $\sigma_0$ back to its proper value, using the same function. As a + // result of this approach, we need to leave sigma_0 as the only + // non-const member variable of this class. + template + class ConstitutiveLaw + { + public: + ConstitutiveLaw (const double E, + const double nu, + const double sigma_0, + const double gamma); + + void + set_sigma_0 (double sigma_zero); + + bool + get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor) const; + + bool + get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, + const std::vector > &point_hessian, + Tensor<5, dim> &stress_strain_tensor_grad) const; + + void + get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor_linearized, + SymmetricTensor<4, dim> &stress_strain_tensor) const; + + private: + const double kappa; + const double mu; + double sigma_0; + const double gamma; + + const SymmetricTensor<4, dim> stress_strain_tensor_kappa; + const SymmetricTensor<4, dim> stress_strain_tensor_mu; + }; + + // The constructor of the ConstitutiveLaw class sets the required material + // parameter for our deformable body. Material parameters for elastic + // isotropic media can be defined in a variety of ways, such as the pair $E, + // \nu$ (elastic modulus and Poisson's number), using the Lame parameters + // $\lambda,mu$ or several other commonly used conventions. Here, the + // constructor takes a description of material parameters in the form of + // $E,\nu$, but since this turns out to these are not the coefficients that + // appear in the equations of the plastic projector, we immediately convert + // them into the more suitable set $\kappa,\mu$ of bulk and shear moduli. In + // addition, the constructor takes $\sigma_0$ (the yield stress absent any + // plastic strain) and $\gamma$ (the hardening parameter) as arguments. In + // this constructor, we also compute the two principal components of the + // stress-strain relation and its linearization. + template + ConstitutiveLaw::ConstitutiveLaw (double E, + double nu, + double sigma_0, + double gamma) + : + //-------------------- + // Plane stress +// kappa (((E*(1+2*nu)) / (std::pow((1+nu),2))) / (3 * (1 - 2 * (nu / (1+nu))))), +// mu (((E*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))), + //-------------------- + // 3d and plane strain + kappa (E / (3 * (1 - 2 * nu))), + mu (E / (2 * (1 + nu))), + //-------------------- + sigma_0(sigma_0), + gamma(gamma), + stress_strain_tensor_kappa (kappa + * outer_product(unit_symmetric_tensor(), + unit_symmetric_tensor())), + stress_strain_tensor_mu (2 * mu + * (identity_tensor() + - outer_product(unit_symmetric_tensor(), + unit_symmetric_tensor()) / 3.0)) + {} + + + template + void + ConstitutiveLaw::set_sigma_0 (double sigma_zero) + { + sigma_0 = sigma_zero; + } + + + // @sect4{ConstitutiveLaw::get_stress_strain_tensor} + + // This is the principal component of the constitutive law. It projects the + // deviatoric part of the stresses in a quadrature point back to the yield + // stress (i.e., the original yield stress $\sigma_0$ plus the term that + // describes linear isotropic hardening). We need this function to calculate + // the nonlinear residual in PlasticityContactProblem::residual_nl_system. The + // computations follow the formulas laid out in the introduction. + // + // The function returns whether the quadrature point is plastic to allow for + // some statistics downstream on how many of the quadrature points are + // plastic and how many are elastic. + template + bool + ConstitutiveLaw:: + get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor) const + { + SymmetricTensor<2, dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) + * strain_tensor; + +// const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor); +// const double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); + const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor); + + stress_strain_tensor = stress_strain_tensor_mu; + if (von_Mises_stress > sigma_0) + { + const double beta = sigma_0 / von_Mises_stress; + stress_strain_tensor *= (gamma + (1 - gamma) * beta); + } + + stress_strain_tensor += stress_strain_tensor_kappa; + + return (von_Mises_stress > sigma_0); + } + + + template + bool + ConstitutiveLaw:: + get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, + const std::vector > &point_hessian, + Tensor<5, dim> &stress_strain_tensor_grad) const + { + SymmetricTensor<2, dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) + * strain_tensor; + + const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor); + const double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); + const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor); + + if (von_Mises_stress > sigma_0) + { + const SymmetricTensor<2, dim> deviator_strain_tensor = deviator(strain_tensor); + const double deviator_strain_tensor_norm = deviator_strain_tensor.norm(); + const double multiplier = -(1-gamma)*sigma_0/(2*mu*std::pow(deviator_strain_tensor_norm,3)); + + Vector multiplier_vector(dim); + multiplier_vector = 0; + + for (unsigned int i=0; i!=dim; ++i) + for (unsigned int m=0; m!=dim; ++m) + for (unsigned int n=0; n!=dim; ++n) + { + multiplier_vector(i) += deviator_strain_tensor[m][n] * + ( 0.5*( point_hessian[m][n][i] + point_hessian[n][m][i] ) + + ( m==n && dim==2 ? -1/dim*(point_hessian[0][0][i] + + point_hessian[1][1][i]) : 0 ) + + ( m==n && dim==3 ? -1/dim*(point_hessian[0][0][i] + + point_hessian[1][1][i] + + point_hessian[2][2][i]) : 0 ) ); + } + + // ----------------------------------------------- + // "Perforated_strip_tension" + // plane stress +// const double VM_factor = std::sqrt(2); + // ----------------------------------------------- + // otherwise + // plane strain / 3d case + const double VM_factor = std::sqrt(1.5); + // ----------------------------------------------- + + for (unsigned int i=0; i!=dim; ++i) + for (unsigned int j=0; j!=dim; ++j) + for (unsigned int k=0; k!=dim; ++k) + for (unsigned int l=0; l!=dim; ++l) + for (unsigned int m=0; m!=dim; ++m) + { + stress_strain_tensor_grad[i][j][k][l][m] = 1/VM_factor + * multiplier + * stress_strain_tensor_mu[i][j][k][l] + * multiplier_vector(m); + } + + }else + { + stress_strain_tensor_grad = 0; + } + + return (von_Mises_stress > sigma_0); + } + + + // @sect4{ConstitutiveLaw::get_linearized_stress_strain_tensors} + + // This function returns the linearized stress strain tensor, linearized + // around the solution $u^{i-1}$ of the previous Newton step $i-1$. The + // parameter strain_tensor (commonly denoted + // $\varepsilon(u^{i-1})$) must be passed as an argument, and serves as the + // linearization point. The function returns the derivative of the nonlinear + // constitutive law in the variable stress_strain_tensor, as well as the + // stress-strain tensor of the linearized problem in + // stress_strain_tensor_linearized. See + // PlasticityContactProblem::assemble_nl_system where this function is used. + template + void + ConstitutiveLaw:: + get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor_linearized, + SymmetricTensor<4, dim> &stress_strain_tensor) const + { + SymmetricTensor<2, dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) + * strain_tensor; + + stress_strain_tensor = stress_strain_tensor_mu; + stress_strain_tensor_linearized = stress_strain_tensor_mu; + + SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor); + const double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); + const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor); + + if (von_Mises_stress > sigma_0) + { + const double beta = sigma_0 / von_Mises_stress; + stress_strain_tensor *= (gamma + (1 - gamma) * beta); + stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta); + deviator_stress_tensor /= deviator_stress_tensor_norm; + stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu + * outer_product(deviator_stress_tensor, + deviator_stress_tensor); + } + + stress_strain_tensor += stress_strain_tensor_kappa; + stress_strain_tensor_linearized += stress_strain_tensor_kappa; + } + + // Finally, below we will need a function that computes the rotation matrix + // induced by a displacement at a given point. In fact, of course, the + // displacement at a single point only has a direction and a magnitude, it + // is the change in direction and magnitude that induces rotations. In + // effect, the rotation matrix can be computed from the gradients of a + // displacement, or, more specifically, from the curl. + // + // The formulas by which the rotation matrices are determined are a little + // awkward, especially in 3d. For 2d, there is a simpler way, so we + // implement this function twice, once for 2d and once for 3d, so that we + // can compile and use the program in both space dimensions if so desired -- + // after all, deal.II is all about dimension independent programming and + // reuse of algorithm thoroughly tested with cheap computations in 2d, for + // the more expensive computations in 3d. Here is one case, where we have to + // implement different algorithms for 2d and 3d, but then can write the rest + // of the program in a way that is independent of the space dimension. + // + // So, without further ado to the 2d implementation: + Tensor<2,2> + get_rotation_matrix (const std::vector > &grad_u) + { + // First, compute the curl of the velocity field from the gradients. Note + // that we are in 2d, so the rotation is a scalar: + const double curl = (grad_u[1][0] - grad_u[0][1]); + + // From this, compute the angle of rotation: + const double angle = std::atan (curl); + + // And from this, build the antisymmetric rotation matrix: + const double t[2][2] = {{ cos(angle), sin(angle) }, + {-sin(angle), cos(angle) } + }; + return Tensor<2,2>(t); + } + + + // The 3d case is a little more contrived: + Tensor<2,3> + get_rotation_matrix (const std::vector > &grad_u) + { + // Again first compute the curl of the velocity field. This time, it is a + // real vector: + const Point<3> curl (grad_u[2][1] - grad_u[1][2], + grad_u[0][2] - grad_u[2][0], + grad_u[1][0] - grad_u[0][1]); + + // From this vector, using its magnitude, compute the tangent of the angle + // of rotation, and from it the actual angle: + const double tan_angle = std::sqrt(curl*curl); + const double angle = std::atan (tan_angle); + + // Now, here's one problem: if the angle of rotation is too small, that + // means that there is no rotation going on (for example a translational + // motion). In that case, the rotation matrix is the identity matrix. + // + // The reason why we stress that is that in this case we have that + // tan_angle==0. Further down, we need to divide by that + // number in the computation of the axis of rotation, and we would get + // into trouble when dividing doing so. Therefore, let's shortcut this and + // simply return the identity matrix if the angle of rotation is really + // small: + if (angle < 1e-9) + { + static const double rotation[3][3] + = {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } }; + static const Tensor<2,3> rot(rotation); + return rot; + } + + // Otherwise compute the real rotation matrix. The algorithm for this is + // not exactly obvious, but can be found in a number of books, + // particularly on computer games where rotation is a very frequent + // operation. Online, you can find a description at + // http://www.makegames.com/3drotation/ and (this particular form, with + // the signs as here) at + // http://www.gamedev.net/reference/articles/article1199.asp: + const double c = std::cos(angle); + const double s = std::sin(angle); + const double t = 1-c; + + const Point<3> axis = curl/tan_angle; + const double rotation[3][3] + = {{ + t *axis[0] *axis[0]+c, + t *axis[0] *axis[1]+s *axis[2], + t *axis[0] *axis[2]-s *axis[1] + }, + { + t *axis[0] *axis[1]-s *axis[2], + t *axis[1] *axis[1]+c, + t *axis[1] *axis[2]+s *axis[0] + }, + { + t *axis[0] *axis[2]+s *axis[1], + t *axis[1] *axis[1]-s *axis[0], + t *axis[2] *axis[2]+c + } + }; + return Tensor<2,3>(rotation); + } + + + //

Equation data: Body forces, boundary forces, + // incremental boundary values

+ // + // The following should be relatively standard. We need classes for + // the boundary forcing term (which we here choose to be zero) + // and incremental boundary values. + namespace EquationData + { + + /* + template + class BoundaryForce : public Function + { + public: + BoundaryForce (); + + virtual + double value (const Point &p, + const unsigned int component = 0) const; + + virtual + void vector_value (const Point &p, + Vector &values) const; + }; + + template + BoundaryForce::BoundaryForce () + : + Function(dim) + {} + + + template + double + BoundaryForce::value (const Point &, + const unsigned int) const + { + return 0.; + } + + template + void + BoundaryForce::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c = 0; c < this->n_components; ++c) + values(c) = BoundaryForce::value(p, c); + } + + // @sect3{The BodyForce class} + // Body forces are generally mediated by one of the four basic + // physical types of forces: + // gravity, strong and weak interaction, and electromagnetism. Unless one + // wants to consider subatomic objects (for which quasistatic deformation is + // irrelevant and an inappropriate description anyway), only gravity and + // electromagnetic forces need to be considered. Let us, for simplicity + // assume that our body has a certain mass density, but is either + // non-magnetic and not electrically conducting or that there are no + // significant electromagnetic fields around. In that case, the body forces + // are simply rho g, where rho is the material + // density and g is a vector in negative z-direction with + // magnitude 9.81 m/s^2. Both the density and g are defined in + // the function, and we take as the density 7700 kg/m^3, a value commonly + // assumed for steel. + // + // To be a little more general and to be able to do computations in 2d as + // well, we realize that the body force is always a function returning a + // dim dimensional vector. We assume that gravity acts along + // the negative direction of the last, i.e. dim-1th + // coordinate. The rest of the implementation of this function should be + // mostly self-explanatory given similar definitions in previous example + // programs. Note that the body force is independent of the location; to + // avoid compiler warnings about unused function arguments, we therefore + // comment out the name of the first argument of the + // vector_value function: + template + class BodyForce : public Function + { + public: + BodyForce (); + + virtual + void + vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + }; + + + template + BodyForce::BodyForce () + : + Function (dim) + {} + + + template + inline + void + BodyForce::vector_value (const Point &p, + Vector &values) const + { + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + + const double g = 9.81; + const double rho = 7700; + + values = 0; + values(dim-1) = -rho * g; + } + + + + template + void + BodyForce::vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + + // @sect3{The IncrementalBoundaryValue class} + + // In addition to body forces, movement can be induced by boundary forces + // and forced boundary displacement. The latter case is equivalent to forces + // being chosen in such a way that they induce certain displacement. + // + // For quasistatic displacement, typical boundary forces would be pressure + // on a body, or tangential friction against another body. We chose a + // somewhat simpler case here: we prescribe a certain movement of (parts of) + // the boundary, or at least of certain components of the displacement + // vector. We describe this by another vector-valued function that, for a + // given point on the boundary, returns the prescribed displacement. + // + // Since we have a time-dependent problem, the displacement increment of the + // boundary equals the displacement accumulated during the length of the + // timestep. The class therefore has to know both the present time and the + // length of the present time step, and can then approximate the incremental + // displacement as the present velocity times the present timestep. + // + // For the purposes of this program, we choose a simple form of boundary + // displacement: we displace the top boundary with constant velocity + // downwards. The rest of the boundary is either going to be fixed (and is + // then described using an object of type ZeroFunction) or free + // (Neumann-type, in which case nothing special has to be done). The + // implementation of the class describing the constant downward motion + // should then be obvious using the knowledge we gained through all the + // previous example programs: + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double present_timestep); + + virtual + void + vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + private: + const double velocity; + const double present_time; + const double present_timestep; + }; + + + template + IncrementalBoundaryValues:: + IncrementalBoundaryValues (const double present_time, + const double present_timestep) + : + Function (dim), + velocity (.1), + present_time (present_time), + present_timestep (present_timestep) + {} + + + template + void + IncrementalBoundaryValues:: + vector_value (const Point &p, + Vector &values) const + { + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + + values = 0; + values(2) = -present_timestep * velocity; + } + + + + template + void + IncrementalBoundaryValues:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + */ + + // ----------------------------- TimoshenkoBeam --------------------------------------- + /* + template + class IncrementalBoundaryForce : public Function + { + public: + IncrementalBoundaryForce (const double present_time, + const double end_time); + + virtual + void vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + private: + const double present_time, + end_time, + shear_force, + length, + depth, + thickness; + }; + + template + IncrementalBoundaryForce:: + IncrementalBoundaryForce (const double present_time, + const double end_time) + : + Function(dim), + present_time (present_time), + end_time (end_time), + shear_force (2e4), + length (.48), + depth (.12), + thickness (.01) + {} + + template + void + IncrementalBoundaryForce::vector_value (const Point &p, + Vector &values) const + { + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + AssertThrow (dim == 2, ExcNotImplemented()); + + // compute traction on the right face of Timoshenko beam problem, t_bar + double inertia_moment = (thickness*std::pow(depth,3)) / 12; + + double x = p(0); + double y = p(1); + + AssertThrow(std::fabs(x-length)<1e-12, ExcNotImplemented()); + + values(0) = 0; + values(1) = - shear_force/(2*inertia_moment) * ( depth*depth/4-y*y ); + + // compute the fraction of imposed force + const double frac = present_time/end_time; + + values *= frac; + } + + template + void + IncrementalBoundaryForce:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + + + template + class BodyForce : public ZeroFunction + { + public: + BodyForce () : ZeroFunction (dim) {} + }; + + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double end_time); + + virtual + void + vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + private: + const double present_time, + end_time, + shear_force, + Youngs_modulus, + Poissons_ratio, + length, + depth, + thickness; + }; + + + template + IncrementalBoundaryValues:: + IncrementalBoundaryValues (const double present_time, + const double end_time) + : + Function (dim), + present_time (present_time), + end_time (end_time), + shear_force (2e4), + Youngs_modulus (2.e11), + Poissons_ratio (.3), + length (.48), + depth (.12), + thickness (.01) + {} + + + template + void + IncrementalBoundaryValues:: + vector_value (const Point &p, + Vector &values) const + { + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + AssertThrow (dim == 2, ExcNotImplemented()); + + + // compute exact displacement of Timoshenko beam problem, u_bar + double inertia_moment = (thickness*std::pow(depth,3)) / 12; + + double x = p(0); + double y = p(1); + + double fac = shear_force / (6*Youngs_modulus*inertia_moment); + + values(0) = fac * y * ( (6*length-3*x)*x + (2+Poissons_ratio)*(y*y-depth*depth/4) ); + values(1) = -fac* ( 3*Poissons_ratio*y*y*(length-x) + 0.25*(4+5*Poissons_ratio)*depth*depth*x + (3*length-x)*x*x ); + + // compute the fraction of imposed force + const double frac = present_time/end_time; + + values *= frac; + } + + + + template + void + IncrementalBoundaryValues:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + */ + + // ------------------------- Thick_tube_internal_pressure ---------------------------------- + /* + template + class IncrementalBoundaryForce : public Function + { + public: + IncrementalBoundaryForce (const double present_time, + const double end_time); + + virtual + void vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + private: + const double present_time, + end_time, + pressure, + inner_radius; + }; + + template + IncrementalBoundaryForce:: + IncrementalBoundaryForce (const double present_time, + const double end_time) + : + Function(dim), + present_time (present_time), + end_time (end_time), + pressure (0.6*2.4e8), +// pressure (1.94e8), + inner_radius(.1) + {} + + template + void + IncrementalBoundaryForce::vector_value (const Point &p, + Vector &values) const + { + AssertThrow (dim == 2, ExcNotImplemented()); + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + + const double eps = 1.e-7 * inner_radius, + radius = p.norm(); + // compute traction on the inner boundary, t_bar + AssertThrow(radius < (eps+inner_radius), ExcInternalError()); + + const double theta = std::atan2(p(1),p(0)); + + values(0) = pressure * std::cos(theta); + values(1) = pressure * std::sin(theta); + + // compute the fraction of imposed force + const double frac = present_time/end_time; + + values *= frac; + } + + template + void + IncrementalBoundaryForce:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + + + template + class BodyForce : public ZeroFunction + { + public: + BodyForce () : ZeroFunction (dim) {} + }; + + + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double end_time); + + virtual + void + vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + private: + const double present_time, + end_time; + }; + + + template + IncrementalBoundaryValues:: + IncrementalBoundaryValues (const double present_time, + const double end_time) + : + Function (dim), + present_time (present_time), + end_time (end_time) + {} + + + template + void + IncrementalBoundaryValues:: + vector_value (const Point &p, + Vector &values) const + { + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + AssertThrow (dim == 2, ExcNotImplemented()); + + values = 0.; + } + + + + template + void + IncrementalBoundaryValues:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + */ + + // ------------------------- Perforated_strip_tension ---------------------------------- + /* + template + class IncrementalBoundaryForce : public Function + { + public: + IncrementalBoundaryForce (const double present_time, + const double end_time); + + virtual + void vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + private: + const double present_time, + end_time; + }; + + template + IncrementalBoundaryForce:: + IncrementalBoundaryForce (const double present_time, + const double end_time) + : + Function(dim), + present_time (present_time), + end_time (end_time) + {} + + template + void + IncrementalBoundaryForce::vector_value (const Point &p, + Vector &values) const + { + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + + values = 0; + + // compute the fraction of imposed force + const double frac = present_time/end_time; + + values *= frac; + } + + template + void + IncrementalBoundaryForce:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + + + template + class BodyForce : public ZeroFunction + { + public: + BodyForce () : ZeroFunction (dim) {} + }; + + + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double end_time); + + virtual + void + vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + private: + const double present_time, + end_time, + imposed_displacement, + height; + }; + + + template + IncrementalBoundaryValues:: + IncrementalBoundaryValues (const double present_time, + const double end_time) + : + Function (dim), + present_time (present_time), + end_time (end_time), + imposed_displacement (0.00055), + height (0.18) + {} + + + template + void + IncrementalBoundaryValues:: + vector_value (const Point &p, + Vector &values) const + { + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + + const double eps = 1.e-8 * height; + + values = 0.; + + // impose displacement only on the top edge + if (std::abs(p[1]-height) < eps) + { + // compute the fraction of imposed displacement + const double inc_frac = 1/end_time; + + values(1) = inc_frac*imposed_displacement; + } + + } + + + + template + void + IncrementalBoundaryValues:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + */ + + // ------------------------- Cantiliver_beam_3d ---------------------------------- + template + class IncrementalBoundaryForce : public Function + { + public: + IncrementalBoundaryForce (const double present_time, + const double end_time); + + virtual + void vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + private: + const double present_time, + end_time, + pressure, + height; + }; + + template + IncrementalBoundaryForce:: + IncrementalBoundaryForce (const double present_time, + const double end_time) + : + Function(dim), + present_time (present_time), + end_time (end_time), + pressure (6e6), + height (200e-3) + {} + + template + void + IncrementalBoundaryForce::vector_value (const Point &p, + Vector &values) const + { + AssertThrow (dim == 3, ExcNotImplemented()); + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + + const double eps = 1.e-7 * height; + + // pressure should be imposed on the top surface, y = height + AssertThrow(std::abs(p[1]-(height/2)) < eps, ExcInternalError()); + + values = 0; + + values(1) = -pressure; + + // compute the fraction of imposed force + const double frac = present_time/end_time; + + values *= frac; + } + + template + void + IncrementalBoundaryForce:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], value_list[p]); + } + + + template + class BodyForce : public ZeroFunction + { + public: + BodyForce () : ZeroFunction (dim) {} + }; + + + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double end_time); + + virtual + void + vector_value (const Point &p, + Vector &values) const; + + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + private: + const double present_time, + end_time; + }; + + + template + IncrementalBoundaryValues:: + IncrementalBoundaryValues (const double present_time, + const double end_time) + : + Function (dim), + present_time (present_time), + end_time (end_time) + {} + + + template + void + IncrementalBoundaryValues:: + vector_value (const Point &p, + Vector &values) const + { + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + AssertThrow (dim == 3, ExcNotImplemented()); + + values = 0.; + } + + + template + void + IncrementalBoundaryValues:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], value_list[p]); + } + + // ------------------------------------------------------------------------------- + } + + + namespace DualFunctional + { + + template + class DualFunctionalBase : public Subscriptor + { + public: + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const = 0; + }; + + + template + class PointValuesEvaluation : public DualFunctionalBase + { + public: + PointValuesEvaluation (const Point &evaluation_point); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + + protected: + const Point evaluation_point; + }; + + + template + PointValuesEvaluation:: + PointValuesEvaluation (const Point &evaluation_point) + : + evaluation_point (evaluation_point) + {} + + + template + void + PointValuesEvaluation:: + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const + { + rhs_dual.reinit (dof_handler_dual.n_dofs()); + const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex; + + typename DoFHandler::active_cell_iterator + cell_dual = dof_handler_dual.begin_active(), + endc_dual = dof_handler_dual.end(); + for (; cell_dual!=endc_dual; ++cell_dual) + for (unsigned int vertex=0; + vertex::vertices_per_cell; + ++vertex) + if (cell_dual->vertex(vertex).distance(evaluation_point) + < cell_dual->diameter()*1e-8) + { + for (unsigned int id=0; id!=dofs_per_vertex; ++id) + { + rhs_dual(cell_dual->vertex_dof_index(vertex,id)) = 1; + } + return; + } + + AssertThrow (false, ExcEvaluationPointNotFound(evaluation_point)); + } + + + template + class PointXDerivativesEvaluation : public DualFunctionalBase + { + public: + PointXDerivativesEvaluation (const Point &evaluation_point); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + + protected: + const Point evaluation_point; + }; + + + template + PointXDerivativesEvaluation:: + PointXDerivativesEvaluation (const Point &evaluation_point) + : + evaluation_point (evaluation_point) + {} + + + template + void + PointXDerivativesEvaluation:: + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const + { + rhs_dual.reinit (dof_handler_dual.n_dofs()); + const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex; + + QGauss quadrature(4); + FEValues fe_values (dof_handler_dual.get_fe(), quadrature, + update_gradients | + update_quadrature_points | + update_JxW_values); + const unsigned int n_q_points = fe_values.n_quadrature_points; + Assert ( n_q_points==quadrature.size() , ExcInternalError() ); + const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell; + + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + double total_volume = 0; + + typename DoFHandler::active_cell_iterator + cell = dof_handler_dual.begin_active(), + endc = dof_handler_dual.end(); + for (; cell!=endc; ++cell) + if (cell->center().distance(evaluation_point) <= + cell->diameter()) + { + fe_values.reinit (cell); + cell_rhs = 0; + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + for (unsigned int i=0; i 0, + ExcEvaluationPointNotFound(evaluation_point)); + + rhs_dual.scale (1./total_volume); + } + + + + template + class MeanDisplacementFace : public DualFunctionalBase + { + public: + MeanDisplacementFace (const unsigned int face_id, + const std::vector comp_mask); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const; + + protected: + const unsigned int face_id; + const std::vector comp_mask; + }; + + + template + MeanDisplacementFace:: + MeanDisplacementFace (const unsigned int face_id, + const std::vector comp_mask ) + : + face_id (face_id), + comp_mask (comp_mask) + { + AssertThrow(comp_mask.size() == dim, + ExcDimensionMismatch (comp_mask.size(), dim) ); + } + + + template + void + MeanDisplacementFace:: + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const + { + AssertThrow (dim >= 2, ExcNotImplemented()); + + rhs_dual.reinit (dof_handler_dual.n_dofs()); + + const QGauss face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1); + FEFaceValues fe_face_values (dof_handler_dual.get_fe(), face_quadrature, + update_values | update_JxW_values); + + const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex; + const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell; + const unsigned int n_face_q_points = face_quadrature.size(); + + AssertThrow(dofs_per_vertex == dim, + ExcDimensionMismatch (dofs_per_vertex, dim) ); + + std::vector comp_vector(dofs_per_vertex); + for (unsigned int i=0; i!=dofs_per_vertex; ++i) + { + if (comp_mask[i]) + { + comp_vector[i] = 1; + } + } + + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + // bound_size : size of the boundary, in 2d is the length + // and in the 3d case, area + double bound_size = 0.; + + typename DoFHandler::active_cell_iterator + cell = dof_handler_dual.begin_active(), + endc = dof_handler_dual.end(); + bool evaluation_face_found = false; + for (; cell!=endc; ++cell) + { + cell_rhs = 0; + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && + cell->face(face)->boundary_indicator() == face_id) + { + if (!evaluation_face_found) + { + evaluation_face_found = true; + } + fe_face_values.reinit (cell, face); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i + class MeanStressFace : public DualFunctionalBase + { + public: + MeanStressFace (const unsigned int face_id, + const std::vector > &comp_stress); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const; + + protected: + const unsigned int face_id; + const std::vector > comp_stress; + }; + + + template + MeanStressFace:: + MeanStressFace (const unsigned int face_id, + const std::vector > &comp_stress ) + : + face_id (face_id), + comp_stress (comp_stress) + { + AssertThrow(comp_stress.size() == dim, + ExcDimensionMismatch (comp_stress.size(), dim) ); + } + + + template + void + MeanStressFace:: + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const + { + AssertThrow (dim >= 2, ExcNotImplemented()); + + rhs_dual.reinit (dof_handler_dual.n_dofs()); + + const QGauss face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1); + + FEFaceValues fe_face_values (dof_handler.get_fe(), face_quadrature, + update_gradients); + FEFaceValues fe_face_values_dual (dof_handler_dual.get_fe(), face_quadrature, + update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell; + const unsigned int n_face_q_points = face_quadrature.size(); + + std::vector > strain_tensor(n_face_q_points); + SymmetricTensor<4, dim> stress_strain_tensor; + + Vector cell_rhs (dofs_per_cell_dual); + + std::vector local_dof_indices (dofs_per_cell_dual); + + // bound_size : size of the boundary, in 2d is the length + // and in the 3d case, area + double bound_size = 0.; + + bool evaluation_face_found = false; + + typename DoFHandler::active_cell_iterator + cell_dual = dof_handler_dual.begin_active(), + endc_dual = dof_handler_dual.end(), + cell = dof_handler.begin_active(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell_dual!=endc_dual; ++cell_dual, ++cell) + { + cell_rhs = 0; + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell_dual->face(face)->at_boundary() + && + cell_dual->face(face)->boundary_indicator() == face_id) + { + if (!evaluation_face_found) + { + evaluation_face_found = true; + } + + fe_face_values.reinit (cell, face); + fe_face_values_dual.reinit (cell_dual, face); + + fe_face_values[displacement].get_function_symmetric_gradients(solution, + strain_tensor); + + for (unsigned int q_point=0; q_point + stress_phi_i = stress_strain_tensor + * fe_face_values_dual[displacement].symmetric_gradient(i, q_point); + + for (unsigned int k=0; k!=dim; ++k) + { + for (unsigned int l=0; l!=dim; ++l) + { + if ( comp_stress[k][l] == 1 ) + { + cell_rhs(i) += stress_phi_i[k][l] + * + fe_face_values_dual.JxW(q_point); + } + + } + } + + } + + } + + } + } + + cell_dual->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i + class MeanStressDomain : public DualFunctionalBase + { + public: + MeanStressDomain (const std::string &base_mesh, + const std::vector > &comp_stress); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const; + + protected: + const std::string base_mesh; + const std::vector > comp_stress; + }; + + + template + MeanStressDomain:: + MeanStressDomain (const std::string &base_mesh, + const std::vector > &comp_stress ) + : + base_mesh (base_mesh), + comp_stress (comp_stress) + { + AssertThrow(comp_stress.size() == dim, + ExcDimensionMismatch (comp_stress.size(), dim) ); + } + + + template + void + MeanStressDomain:: + assemble_rhs (const DoFHandler &dof_handler, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DoFHandler &dof_handler_dual, + Vector &rhs_dual) const + { + AssertThrow (base_mesh == "Cantiliver_beam_3d", ExcNotImplemented()); + AssertThrow (dim == 3, ExcNotImplemented()); + + // Mean stress at the specified domain is of interest. + // The interest domains are located on the bottom and top of the flanges + // close to the clamped face, z = 0 + // top domain: height/2 - thickness_flange <= y <= height/2 + // 0 <= z <= 2 * thickness_flange + // bottom domain: -height/2 <= y <= -height/2 + thickness_flange + // 0 <= z <= 2 * thickness_flange + + const double height = 200e-3, + thickness_flange = 10e-3; + + rhs_dual.reinit (dof_handler_dual.n_dofs()); + + const QGauss quadrature_formula(dof_handler_dual.get_fe().tensor_degree()+1); + + FEValues fe_values (dof_handler.get_fe(), quadrature_formula, + update_gradients); + FEValues fe_values_dual (dof_handler_dual.get_fe(), quadrature_formula, + update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + std::vector > strain_tensor(n_q_points); + SymmetricTensor<4, dim> stress_strain_tensor; + + Vector cell_rhs (dofs_per_cell_dual); + + std::vector local_dof_indices (dofs_per_cell_dual); + + // domain_size : size of the interested domain, in 2d is the area + // and in the 3d case, volume + double domain_size = 0.; + + bool evaluation_domain_found = false; + + typename DoFHandler::active_cell_iterator + cell_dual = dof_handler_dual.begin_active(), + endc_dual = dof_handler_dual.end(), + cell = dof_handler.begin_active(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell_dual!=endc_dual; ++cell_dual, ++cell) + { + const double y = cell->center()[1], + z = cell->center()[2]; + // top domain: height/2 - thickness_flange <= y <= height/2 + // 0 <= z <= 2 * thickness_flange + // bottom domain: -height/2 <= y <= -height/2 + thickness_flange + // 0 <= z <= 2 * thickness_flange + if ( ((z > 0) && (z < 2*thickness_flange)) && + ( ((y > height/2 - thickness_flange) && (y < height/2)) || + ((y > -height/2) && (y < -height/2 + thickness_flange)) ) ) + { + cell_rhs = 0; + + if (!evaluation_domain_found) + { + evaluation_domain_found = true; + } + + fe_values.reinit(cell); + fe_values_dual.reinit(cell_dual); + + fe_values[displacement].get_function_symmetric_gradients(solution, + strain_tensor); + + for (unsigned int q_point=0; q_point + stress_phi_i = stress_strain_tensor + * fe_values_dual[displacement].symmetric_gradient(i, q_point); + + for (unsigned int k=0; k!=dim; ++k) + { + for (unsigned int l=0; l!=dim; ++l) + { + if ( comp_stress[k][l] == 1 ) + { + cell_rhs(i) += stress_phi_i[k][l] + * + fe_values_dual.JxW(q_point); + } + + } + } + + } + + } + + } + + cell_dual->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i + class MeanStrainEnergyFace : public DualFunctionalBase + { + public: + MeanStrainEnergyFace (const unsigned int face_id, + const Function &lambda_function, + const Function &mu_function ); + + void assemble_rhs_nonlinear (const DoFHandler &primal_dof_handler, + const Vector &primal_solution, + const DoFHandler &dof_handler, + Vector &rhs) const; + + protected: + const unsigned int face_id; + const SmartPointer > lambda_function; + const SmartPointer > mu_function; + }; + + + template + MeanStrainEnergyFace:: + MeanStrainEnergyFace (const unsigned int face_id, + const Function &lambda_function, + const Function &mu_function ) + : + face_id (face_id), + lambda_function (&lambda_function), + mu_function (&mu_function) + {} + + + template + void + MeanStrainEnergyFace:: + assemble_rhs_nonlinear (const DoFHandler &primal_dof_handler, + const Vector &primal_solution, + const DoFHandler &dof_handler, + Vector &rhs) const + { + // Assemble right hand side of the dual problem when the quantity of interest is + // a nonlinear functinoal. In this case, the QoI should be linearized which depends + // on the solution of the primal problem. + // The extracter of the linearized QoI functional is the gradient of the the original + // QoI functional with the primal solution values. + + AssertThrow (dim >= 2, ExcNotImplemented()); + + rhs.reinit (dof_handler.n_dofs()); + + const QGauss face_quadrature(dof_handler.get_fe().tensor_degree()+1); + FEFaceValues primal_fe_face_values (primal_dof_handler.get_fe(), face_quadrature, + update_quadrature_points | + update_gradients | update_hessians | + update_JxW_values); + + FEFaceValues fe_face_values (dof_handler.get_fe(), face_quadrature, + update_values); + + const unsigned int dofs_per_vertex = primal_dof_handler.get_fe().dofs_per_vertex; + const unsigned int n_face_q_points = face_quadrature.size(); + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + + AssertThrow(dofs_per_vertex == dim, + ExcDimensionMismatch (dofs_per_vertex, dim) ); + + std::vector< std::vector< Tensor<1,dim> > > primal_solution_gradients; + primal_solution_gradients.resize(n_face_q_points); + + std::vector > > primal_solution_hessians; + primal_solution_hessians.resize (n_face_q_points); + + for (unsigned int i=0; i!=n_face_q_points; ++i) + { + primal_solution_gradients[i].resize (dofs_per_vertex); + primal_solution_hessians[i].resize (dofs_per_vertex); + } + + std::vector lambda_values (n_face_q_points); + std::vector mu_values (n_face_q_points); + + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + // bound_size : size of the boundary, in 2d is the length + // and in the 3d case, area + double bound_size = 0.; + + bool evaluation_face_found = false; + + typename DoFHandler::active_cell_iterator + primal_cell = primal_dof_handler.begin_active(), + primal_endc = primal_dof_handler.end(); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell, ++primal_cell) + { + cell_rhs = 0; + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && + cell->face(face)->boundary_indicator() == face_id) + { + if (!evaluation_face_found) + { + evaluation_face_found = true; + } + primal_fe_face_values.reinit (primal_cell, face); + + primal_fe_face_values.get_function_grads (primal_solution, + primal_solution_gradients); + + primal_fe_face_values.get_function_hessians (primal_solution, + primal_solution_hessians); + + lambda_function->value_list (primal_fe_face_values.get_quadrature_points(), lambda_values); + mu_function->value_list (primal_fe_face_values.get_quadrature_points(), mu_values); + + fe_face_values.reinit (cell, face); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i + class DualSolver + { + public: + DualSolver (const Triangulation &triangulation, + const FESystem &fe, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DualFunctional::DualFunctionalBase &dual_functional, + const unsigned int ×tep_no, + const std::string &output_dir, + const std::string &base_mesh, + const double &present_time, + const double &end_time); + + void compute_error_DWR (Vector &estimated_error_per_cell); + + ~DualSolver (); + + private: + void setup_system (); + void compute_dirichlet_constraints (); + void assemble_matrix (); + void assemble_rhs (); + void solve (); + void output_results (); + + const FESystem fe; + DoFHandler dof_handler; + const Vector solution; + + const unsigned int fe_degree; + + + const unsigned int fe_degree_dual; + FESystem fe_dual; + DoFHandler dof_handler_dual; + + const QGauss quadrature_formula; + const QGauss face_quadrature_formula; + + ConstraintMatrix constraints_hanging_nodes_dual; + ConstraintMatrix constraints_dirichlet_and_hanging_nodes_dual; + + SparsityPattern sparsity_pattern_dual; + SparseMatrix system_matrix_dual; + Vector system_rhs_dual; + Vector solution_dual; + + const ConstitutiveLaw constitutive_law; + + const SmartPointer > triangulation; + const SmartPointer > dual_functional; + + unsigned int timestep_no; + std::string output_dir; + const std::string base_mesh; + double present_time; + double end_time; + }; + + + template + DualSolver:: + DualSolver (const Triangulation &triangulation, + const FESystem &fe, + const Vector &solution, + const ConstitutiveLaw &constitutive_law, + const DualFunctional::DualFunctionalBase &dual_functional, + const unsigned int ×tep_no, + const std::string &output_dir, + const std::string &base_mesh, + const double &present_time, + const double &end_time) + : + fe (fe), + dof_handler (triangulation), + solution(solution), + fe_degree(fe.tensor_degree()), + fe_degree_dual(fe_degree + 1), + fe_dual(FE_Q(fe_degree_dual), dim), + dof_handler_dual (triangulation), + quadrature_formula (fe_degree_dual + 1), + face_quadrature_formula (fe_degree_dual + 1), + constitutive_law (constitutive_law), + triangulation (&triangulation), + dual_functional (&dual_functional), + timestep_no (timestep_no), + output_dir (output_dir), + base_mesh (base_mesh), + present_time (present_time), + end_time (end_time) + {} + + + template + DualSolver::~DualSolver() + { + dof_handler_dual.clear (); + } + + + template + void DualSolver::setup_system() + { + dof_handler.distribute_dofs(fe); + + dof_handler_dual.distribute_dofs (fe_dual); + std::cout << " Number of degrees of freedom in dual problem: " + << dof_handler_dual.n_dofs() + << std::endl; + + constraints_hanging_nodes_dual.clear (); + DoFTools::make_hanging_node_constraints (dof_handler_dual, + constraints_hanging_nodes_dual); + constraints_hanging_nodes_dual.close (); + + compute_dirichlet_constraints(); + + sparsity_pattern_dual.reinit (dof_handler_dual.n_dofs(), + dof_handler_dual.n_dofs(), + dof_handler_dual.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler_dual, sparsity_pattern_dual); + +// constraints_hanging_nodes_dual.condense (sparsity_pattern_dual); + constraints_dirichlet_and_hanging_nodes_dual.condense (sparsity_pattern_dual); + + sparsity_pattern_dual.compress(); + + system_matrix_dual.reinit (sparsity_pattern_dual); + + solution_dual.reinit (dof_handler_dual.n_dofs()); + system_rhs_dual.reinit (dof_handler_dual.n_dofs()); + + } + + template + void DualSolver::compute_dirichlet_constraints() + { + constraints_dirichlet_and_hanging_nodes_dual.clear (); + constraints_dirichlet_and_hanging_nodes_dual.merge(constraints_hanging_nodes_dual); + + std::vector component_mask(dim); + + if (base_mesh == "Timoshenko beam") + { + VectorTools::interpolate_boundary_values(dof_handler_dual, + 0, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes_dual, + ComponentMask()); + }else if (base_mesh == "Thick_tube_internal_pressure") + { + // the boundary x = 0 + component_mask[0] = true; component_mask[1] = false; + VectorTools::interpolate_boundary_values (dof_handler_dual, + 2, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes_dual, + component_mask); + // the boundary y = 0 + component_mask[0] = false; component_mask[1] = true; + VectorTools::interpolate_boundary_values (dof_handler_dual, + 3, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes_dual, + component_mask); + }else if (base_mesh == "Perforated_strip_tension") + { + // the boundary x = 0 + component_mask[0] = true; component_mask[1] = false; component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler_dual, + 4, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes_dual, + component_mask); + // the boundary y = 0 + component_mask[0] = false; component_mask[1] = true; component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler_dual, + 1, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes_dual, + component_mask); + // the boundary y = imposed incremental displacement + component_mask[0] = false; component_mask[1] = true; component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler_dual, + 3, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes_dual, + component_mask); + }else if (base_mesh == "Cantiliver_beam_3d") + { + // the boundary x = y = z = 0 + component_mask[0] = true; component_mask[1] = true; component_mask[2] = true; + VectorTools::interpolate_boundary_values (dof_handler_dual, + 1, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes_dual, + component_mask); + }else + { + AssertThrow(false, ExcNotImplemented()); + } + + constraints_dirichlet_and_hanging_nodes_dual.close(); + } + + + template + void DualSolver::assemble_matrix() + { + FEValues fe_values(fe, quadrature_formula, update_gradients); + + FEValues fe_values_dual(fe_dual, quadrature_formula, + update_values | update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell_dual = fe_dual.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell_dual, dofs_per_cell_dual); + + std::vector local_dof_indices(dofs_per_cell_dual); + + typename DoFHandler::active_cell_iterator + cell_dual = dof_handler_dual.begin_active(), + endc_dual = dof_handler_dual.end(), + cell = dof_handler.begin_active(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell_dual != endc_dual; ++cell_dual, ++cell) + if (cell_dual->is_locally_owned()) + { + fe_values.reinit(cell); + + fe_values_dual.reinit(cell_dual); + cell_matrix = 0; + + std::vector > strain_tensor(n_q_points); + fe_values[displacement].get_function_symmetric_gradients(solution, + strain_tensor); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + SymmetricTensor<4, dim> stress_strain_tensor_linearized; + SymmetricTensor<4, dim> stress_strain_tensor; + constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point], + stress_strain_tensor_linearized, + stress_strain_tensor); + + for (unsigned int i = 0; i < dofs_per_cell_dual; ++i) + { + const SymmetricTensor<2, dim> + stress_phi_i = stress_strain_tensor_linearized + * fe_values_dual[displacement].symmetric_gradient(i, q_point); + + for (unsigned int j = 0; j < dofs_per_cell_dual; ++j) + cell_matrix(i, j) += (stress_phi_i + * fe_values_dual[displacement].symmetric_gradient(j, q_point) + * fe_values_dual.JxW(q_point)); + + } + + } + + cell_dual->get_dof_indices(local_dof_indices); + constraints_dirichlet_and_hanging_nodes_dual.distribute_local_to_global(cell_matrix, + local_dof_indices, + system_matrix_dual); + + } + + } + + + template + void DualSolver::assemble_rhs() + { + dual_functional->assemble_rhs (dof_handler, solution, constitutive_law, + dof_handler_dual, system_rhs_dual); + constraints_dirichlet_and_hanging_nodes_dual.condense (system_rhs_dual); + } + + + template + void DualSolver::solve() + { + // +++ direct solver +++++++++ + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix_dual); + + // After the decomposition, we can use A_direct like a matrix representing + // the inverse of our system matrix, so to compute the solution we just + // have to multiply with the right hand side vector: + A_direct.vmult(solution_dual, system_rhs_dual); + + // ++++ iterative solver ++ CG ++++ doesn't work +// SolverControl solver_control (5000, 1e-12); +// SolverCG<> cg (solver_control); +// +// PreconditionSSOR<> preconditioner; +// preconditioner.initialize(system_matrix_dual, 1.2); +// +// cg.solve (system_matrix_dual, solution_dual, system_rhs_dual, +// preconditioner); + + // ++++ iterative solver ++ BiCGStab ++++++ doesn't work +// SolverControl solver_control (5000, 1e-12); +// SolverBicgstab<> bicgstab (solver_control); +// +// PreconditionJacobi<> preconditioner; +// preconditioner.initialize(system_matrix_dual, 1.0); +// +// bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual, +// preconditioner); + + // +++++++++++++++++++++++++++++++++++++++++++++++++ + + constraints_dirichlet_and_hanging_nodes_dual.distribute (solution_dual); + } + + template + void DualSolver::output_results() + { + std::string filename = (output_dir + "dual-solution-" + + Utilities::int_to_string(timestep_no, 4) + ".vtk"); + std::ofstream output (filename.c_str()); + DataOut data_out; + data_out.attach_dof_handler (dof_handler_dual); + std::vector solution_names; + switch (dim) + { + case 1: + solution_names.push_back ("displacement"); + break; + case 2: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + break; + case 3: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + solution_names.push_back ("z_displacement"); + break; + default: + Assert (false, ExcNotImplemented()); + } + data_out.add_data_vector (solution_dual, solution_names); + data_out.build_patches (); + data_out.write_vtk (output); + } + + template + void DualSolver::compute_error_DWR (Vector &estimated_error_per_cell) + { + Assert (estimated_error_per_cell.size() == triangulation->n_global_active_cells(), + ExcDimensionMismatch (estimated_error_per_cell.size(), triangulation->n_global_active_cells())); + + // solve the dual problem + setup_system (); + assemble_matrix (); + assemble_rhs (); + solve (); + output_results (); + + // compuate the dual weights + Vector primal_solution (dof_handler_dual.n_dofs()); + FETools::interpolate (dof_handler, + solution, + dof_handler_dual, + constraints_dirichlet_and_hanging_nodes_dual, + primal_solution); + + ConstraintMatrix constraints_hanging_nodes; + DoFTools::make_hanging_node_constraints (dof_handler, + constraints_hanging_nodes); + constraints_hanging_nodes.close(); + Vector dual_weights (dof_handler_dual.n_dofs()); + FETools::interpolation_difference (dof_handler_dual, + constraints_dirichlet_and_hanging_nodes_dual, + solution_dual, + dof_handler, + constraints_hanging_nodes, + dual_weights); + + // estimate the error + FEValues fe_values(fe_dual, quadrature_formula, + update_values | + update_gradients | + update_hessians | + update_quadrature_points | + update_JxW_values); + + const unsigned int n_q_points = quadrature_formula.size(); + std::vector > strain_tensor(n_q_points); + SymmetricTensor<4, dim> stress_strain_tensor_linearized; + SymmetricTensor<4, dim> stress_strain_tensor; + Tensor<5, dim> stress_strain_tensor_grad; + std::vector > > cell_hessians (n_q_points); + for (unsigned int i=0; i!=n_q_points; ++i) + { + cell_hessians[i].resize (dim); + } + std::vector > dual_weights_cell_values (n_q_points, Vector(dim)); + + const EquationData::BodyForce body_force; + std::vector > body_force_values (n_q_points, Vector(dim)); + const FEValuesExtractors::Vector displacement(0); + + + FEFaceValues fe_face_values_cell(fe_dual, face_quadrature_formula, + update_values | + update_quadrature_points| + update_gradients | + update_JxW_values | + update_normal_vectors), + fe_face_values_neighbor (fe_dual, face_quadrature_formula, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors); + FESubfaceValues fe_subface_values_cell (fe_dual, face_quadrature_formula, + update_gradients); + + const unsigned int n_face_q_points = face_quadrature_formula.size(); + std::vector > jump_residual (n_face_q_points, Vector(dim)); + std::vector > dual_weights_face_values (n_face_q_points, Vector(dim)); + + std::vector > > cell_grads(n_face_q_points); + for (unsigned int i=0; i!=n_face_q_points; ++i) + { + cell_grads[i].resize (dim); + } + std::vector > > neighbor_grads(n_face_q_points); + for (unsigned int i=0; i!=n_face_q_points; ++i) + { + neighbor_grads[i].resize (dim); + } + SymmetricTensor<2, dim> q_cell_strain_tensor; + SymmetricTensor<2, dim> q_neighbor_strain_tensor; + SymmetricTensor<4, dim> cell_stress_strain_tensor; + SymmetricTensor<4, dim> neighbor_stress_strain_tensor; + + + typename std::map::face_iterator, Vector > + face_integrals; + typename DoFHandler::active_cell_iterator + cell = dof_handler_dual.begin_active(), + endc = dof_handler_dual.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + { + face_integrals[cell->face(face_no)].reinit (dim); + face_integrals[cell->face(face_no)] = -1e20; + } + } + + std::vector > error_indicators_vector; + error_indicators_vector.resize( triangulation->n_active_cells(), + Vector(dim) ); + + // ----------------- estimate_some ------------------------- + cell = dof_handler_dual.begin_active(); + unsigned int present_cell = 0; + for (; cell!=endc; ++cell, ++present_cell) + if (cell->is_locally_owned()) + { + // --------------- integrate_over_cell ------------------- + fe_values.reinit(cell); + body_force.vector_value_list(fe_values.get_quadrature_points(), + body_force_values); + fe_values[displacement].get_function_symmetric_gradients(primal_solution, + strain_tensor); + fe_values.get_function_hessians(primal_solution, cell_hessians); + + fe_values.get_function_values(dual_weights, + dual_weights_cell_values); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point], + stress_strain_tensor_linearized, + stress_strain_tensor); + constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point], + cell_hessians[q_point], + stress_strain_tensor_grad); + + for (unsigned int i=0; i!=dim; ++i) + { + error_indicators_vector[present_cell](i) += + body_force_values[q_point](i)* + dual_weights_cell_values[q_point](i)* + fe_values.JxW(q_point); + for (unsigned int j=0; j!=dim; ++j) + { + for (unsigned int k=0; k!=dim; ++k) + { + for (unsigned int l=0; l!=dim; ++l) + { + error_indicators_vector[present_cell](i) += + ( stress_strain_tensor[i][j][k][l]* + 0.5*(cell_hessians[q_point][k][l][j] + + + cell_hessians[q_point][l][k][j]) + + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l] + ) * + dual_weights_cell_values[q_point](i) * + fe_values.JxW(q_point); + } + } + } + + } + + } + // ------------------------------------------------------- + // compute face_integrals + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + { + if (cell->face(face_no)->at_boundary()) + { + for (unsigned int id=0; id!=dim; ++id) + { + face_integrals[cell->face(face_no)](id) = 0; + } + continue; + } + + if ((cell->neighbor(face_no)->has_children() == false) && + (cell->neighbor(face_no)->level() == cell->level()) && + (cell->neighbor(face_no)->index() < cell->index())) + continue; + + if (cell->at_boundary(face_no) == false) + if (cell->neighbor(face_no)->level() < cell->level()) + continue; + + + if (cell->face(face_no)->has_children() == false) + { + // ------------- integrate_over_regular_face ----------- + fe_face_values_cell.reinit(cell, face_no); + fe_face_values_cell.get_function_grads (primal_solution, + cell_grads); + + Assert (cell->neighbor(face_no).state() == IteratorState::valid, + ExcInternalError()); + const unsigned int + neighbor_neighbor = cell->neighbor_of_neighbor (face_no); + const typename DoFHandler::active_cell_iterator + neighbor = cell->neighbor(face_no); + + fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor); + fe_face_values_neighbor.get_function_grads (primal_solution, + neighbor_grads); + + for (unsigned int q_point=0; q_point face_integral_vector(dim); + face_integral_vector = 0; + for (unsigned int q_point=0; q_pointface(face_no)) != face_integrals.end(), + ExcInternalError()); + + for (unsigned int i=0; i!=dim; ++i) + { + Assert (face_integrals[cell->face(face_no)](i) == -1e20, + ExcInternalError()); + face_integrals[cell->face(face_no)](i) = face_integral_vector(i); + + } + + // ----------------------------------------------------- + }else + { + // ------------- integrate_over_irregular_face --------- + const typename DoFHandler::face_iterator + face = cell->face(face_no); + const typename DoFHandler::cell_iterator + neighbor = cell->neighbor(face_no); + Assert (neighbor.state() == IteratorState::valid, + ExcInternalError()); + Assert (neighbor->has_children(), + ExcInternalError()); + + const unsigned int + neighbor_neighbor = cell->neighbor_of_neighbor (face_no); + + for (unsigned int subface_no=0; + subface_non_children(); ++subface_no) + { + const typename DoFHandler::active_cell_iterator + neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no); + Assert (neighbor_child->face(neighbor_neighbor) == + cell->face(face_no)->child(subface_no), + ExcInternalError()); + + fe_subface_values_cell.reinit (cell, face_no, subface_no); + fe_subface_values_cell.get_function_grads (primal_solution, + cell_grads); + fe_face_values_neighbor.reinit (neighbor_child, + neighbor_neighbor); + fe_face_values_neighbor.get_function_grads (primal_solution, + neighbor_grads); + + for (unsigned int q_point=0; q_point face_integral_vector(dim); + face_integral_vector = 0; + for (unsigned int q_point=0; q_pointface(neighbor_neighbor)](i) = face_integral_vector(i); + } + + } + + Vector sum (dim); + sum = 0; + for (unsigned int subface_no=0; + subface_non_children(); ++subface_no) + { + Assert (face_integrals.find(face->child(subface_no)) != + face_integrals.end(), + ExcInternalError()); + for (unsigned int i=0; i!=dim; ++i) + { + Assert (face_integrals[face->child(subface_no)](i) != -1e20, + ExcInternalError()); + sum(i) += face_integrals[face->child(subface_no)](i); + } + } + for (unsigned int i=0; i!=dim; ++i) + { + face_integrals[face](i) = sum(i); + } + + + // ----------------------------------------------------- + } + + + } + } + // ---------------------------------------------------------- + + present_cell=0; + cell = dof_handler_dual.begin_active(); + for (; cell!=endc; ++cell, ++present_cell) + if (cell->is_locally_owned()) + { + for (unsigned int face_no=0; face_no::faces_per_cell; + ++face_no) + { + Assert(face_integrals.find(cell->face(face_no)) != + face_integrals.end(), + ExcInternalError()); + + for (unsigned int id=0; id!=dim; ++id) + { + error_indicators_vector[present_cell](id) + -= 0.5*face_integrals[cell->face(face_no)](id); + } + + } + + estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm(); + + } + } + + + + // @sect3{The PlasticityContactProblem class template} + + // This is the main class of this program and supplies all functions + // and variables needed to describe + // the nonlinear contact problem. It is + // close to step-41 but with some additional + // features like handling hanging nodes, + // a Newton method, using Trilinos and p4est + // for parallel distributed computing. + // To deal with hanging nodes makes + // life a bit more complicated since + // we need another ConstraintMatrix now. + // We create a Newton method for the + // active set method for the contact + // situation and to handle the nonlinear + // operator for the constitutive law. + // + // The general layout of this class is very much like for most other tutorial programs. + // To make our life a bit easier, this class reads a set of input parameters from an input file. These + // parameters, using the ParameterHandler class, are declared in the declare_parameters + // function (which is static so that it can be called before we even create an object of the current + // type), and a ParameterHandler object that has been used to read an input file will then be passed + // to the constructor of this class. + // + // The remaining member functions are by and large as we have seen in several of the other tutorial + // programs, though with additions for the current nonlinear system. We will comment on their purpose + // as we get to them further below. + template + class ElastoPlasticProblem + { + public: + ElastoPlasticProblem (const ParameterHandler &prm); + + void run (); + + static void declare_parameters (ParameterHandler &prm); + + private: + void make_grid (); + void setup_system (); + void compute_dirichlet_constraints (); + void assemble_newton_system (const TrilinosWrappers::MPI::Vector &linearization_point, + const TrilinosWrappers::MPI::Vector &delta_linearization_point); + void compute_nonlinear_residual (const TrilinosWrappers::MPI::Vector &linearization_point); + void solve_newton_system (); + void solve_newton (); + void compute_error (); + void compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution); + void refine_grid (); + void move_mesh (const TrilinosWrappers::MPI::Vector &displacement) const; + void output_results (const std::string &filename_base); + + // Next are three functions that handle the history variables stored in each + // quadrature point. The first one is called before the first timestep to + // set up a pristine state for the history variables. It only works on + // those quadrature points on cells that belong to the present processor: + void setup_quadrature_point_history (); + + // The second one updates the history variables at the end of each + // timestep: + void update_quadrature_point_history (); + + // As far as member variables are concerned, we start with ones that we use to + // indicate the MPI universe this program runs on, and then two numbers + // telling us how many participating processors there are, and where in + // this world we are., a stream we use to let + // exactly one processor produce output to the console (see step-17) and + // a variable that is used to time the various sections of the program: + MPI_Comm mpi_communicator; + const unsigned int n_mpi_processes; + const unsigned int this_mpi_process; + ConditionalOStream pcout; + TimerOutput computing_timer; + + // The next group describes the mesh and the finite element space. + // In particular, for this parallel program, the finite element + // space has associated with it variables that indicate which degrees + // of freedom live on the current processor (the index sets, see + // also step-40 and the @ref distributed documentation module) as + // well as a variety of constraints: those imposed by hanging nodes, + // by Dirichlet boundary conditions, and by the active set of + // contact nodes. Of the three ConstraintMatrix variables defined + // here, the first only contains hanging node constraints, the + // second also those associated with Dirichlet boundary conditions, + // and the third these plus the contact constraints. + // + // The variable active_set consists of those degrees + // of freedom constrained by the contact, and we use + // fraction_of_plastic_q_points_per_cell to keep + // track of the fraction of quadrature points on each cell where + // the stress equals the yield stress. The latter is only used to + // create graphical output showing the plastic zone, but not for + // any further computation; the variable is a member variable of + // this class since the information is computed as a by-product + // of computing the residual, but is used only much later. (Note + // that the vector is a vector of length equal to the number of + // active cells on the local mesh; it is never used to + // exchange information between processors and can therefore be + // a regular deal.II vector.) + const unsigned int n_initial_global_refinements; + parallel::distributed::Triangulation triangulation; + + const unsigned int fe_degree; + FESystem fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + ConstraintMatrix constraints_hanging_nodes; + ConstraintMatrix constraints_dirichlet_and_hanging_nodes; + + Vector fraction_of_plastic_q_points_per_cell; + + // One difference of this program is that we declare the quadrature + // formula in the class declaration. The reason is that in all the other + // programs, it didn't do much harm if we had used different quadrature + // formulas when computing the matrix and the right hand side, for + // example. However, in the present case it does: we store information in + // the quadrature points, so we have to make sure all parts of the program + // agree on where they are and how many there are on each cell. Thus, let + // us first declare the quadrature formula that will be used throughout... + const QGauss quadrature_formula; + const QGauss face_quadrature_formula; + + // ... and then also have a vector of history objects, one per quadrature + // point on those cells for which we are responsible (i.e. we don't store + // history data for quadrature points on cells that are owned by other + // processors). + std::vector > quadrature_point_history; + + // The way this object is accessed is through a user pointer + // that each cell, face, or edge holds: it is a void* pointer + // that can be used by application programs to associate arbitrary data to + // cells, faces, or edges. What the program actually does with this data + // is within its own responsibility, the library just allocates some space + // for these pointers, and application programs can set and read the + // pointers for each of these objects. + + + // The next block of variables corresponds to the solution + // and the linear systems we need to form. In particular, this + // includes the Newton matrix and right hand side; the vector + // that corresponds to the residual (i.e., the Newton right hand + // side) but from which we have not eliminated the various + // constraints and that is used to determine which degrees of + // freedom need to be constrained in the next iteration; and + // a vector that corresponds to the diagonal of the $B$ matrix + // briefly mentioned in the introduction and discussed in the + // accompanying paper. + TrilinosWrappers::SparseMatrix newton_matrix; + + TrilinosWrappers::MPI::Vector solution; + TrilinosWrappers::MPI::Vector incremental_displacement; + TrilinosWrappers::MPI::Vector newton_rhs; + TrilinosWrappers::MPI::Vector newton_rhs_residual; + + // The next block of variables is then related to the time dependent + // nature of the problem: they denote the length of the time interval + // which we want to simulate, the present time and number of time step, + // and length of present timestep: + double present_time; + double present_timestep; + double end_time; + unsigned int timestep_no; + + // The next block contains the variables that describe the material + // response: + const double e_modulus, nu, sigma_0, gamma; + ConstitutiveLaw constitutive_law; + + // And then there is an assortment of other variables that are used + // to identify the mesh we are asked to build as selected by the + // parameter file, the obstacle that is being pushed into the + // deformable body, the mesh refinement strategy, whether to transfer + // the solution from one mesh to the next, and how many mesh + // refinement cycles to perform. As possible, we mark these kinds + // of variables as const to help the reader identify + // which ones may or may not be modified later on (the output directory + // being an exception -- it is never modified outside the constructor + // but it is awkward to initialize in the member-initializer-list + // following the colon in the constructor since there we have only + // one shot at setting it; the same is true for the mesh refinement + // criterion): + const std::string base_mesh; + + struct RefinementStrategy + { + enum value + { + refine_global, + refine_percentage, + refine_fix_dofs + }; + }; + typename RefinementStrategy::value refinement_strategy; + + struct ErrorEstimationStrategy + { + enum value + { + kelly_error, + residual_error, + weighted_residual_error, + weighted_kelly_error + }; + }; + typename ErrorEstimationStrategy::value error_estimation_strategy; + + Vector estimated_error_per_cell; + + const bool transfer_solution; + std::string output_dir; + TableHandler table_results, + table_results_2, + table_results_3; + + unsigned int current_refinement_cycle; + + const double max_relative_error; + float relative_error; + + const bool show_stresses; + }; + + + // @sect3{Implementation of the PlasticityContactProblem class} + + // @sect4{PlasticityContactProblem::declare_parameters} + + // Let us start with the declaration of run-time parameters that can be + // selected in the input file. These values will be read back in the + // constructor of this class to initialize the member variables of this + // class: + template + void + ElastoPlasticProblem::declare_parameters (ParameterHandler &prm) + { + prm.declare_entry("polynomial degree", "1", + Patterns::Integer(), + "Polynomial degree of the FE_Q finite element space, typically 1 or 2."); + prm.declare_entry("number of initial refinements", "2", + Patterns::Integer(), + "Number of initial global mesh refinement steps before " + "the first computation."); + prm.declare_entry("refinement strategy", "percentage", + Patterns::Selection("global|percentage"), + "Mesh refinement strategy:\n" + " global: one global refinement\n" + " percentage: a fixed percentage of cells gets refined using the selected error estimator."); + prm.declare_entry("error estimation strategy", "kelly_error", + Patterns::Selection("kelly_error|residual_error|weighted_residual_error"), + "Error estimation strategy:\n" + " kelly_error: Kelly error estimator\n" + " residual_error: residual-based error estimator\n" + " weighted_residual_error: dual weighted residual (Goal-oriented) error estimator.\n"); + prm.declare_entry("maximum relative error","0.05", + Patterns::Double(), + "maximum relative error which plays the role of a criteria for refinement."); + prm.declare_entry("number of cycles", "5", + Patterns::Integer(), + "Number of adaptive mesh refinement cycles to run."); + prm.declare_entry("output directory", "", + Patterns::Anything(), + "Directory for output files (graphical output and benchmark " + "statistics). If empty, use the current directory."); + prm.declare_entry("transfer solution", "true", + Patterns::Bool(), + "Whether the solution should be used as a starting guess " + "for the next finer mesh. If false, then the iteration starts at " + "zero on every mesh."); + prm.declare_entry("base mesh", "Thick_tube_internal_pressure", + Patterns::Selection("Timoshenko beam|Thick_tube_internal_pressure|" + "Perforated_strip_tension|Cantiliver_beam_3d"), + "Select the shape of the domain: 'box' or 'half sphere'"); + prm.declare_entry("elasticity modulus","2.e11", + Patterns::Double(), + "Elasticity modulus of the material in MPa (N/mm2)"); + prm.declare_entry("Poissons ratio","0.3", + Patterns::Double(), + "Poisson's ratio of the material"); + prm.declare_entry("yield stress","2.e11", + Patterns::Double(), + "Yield stress of the material in MPa (N/mm2)"); + prm.declare_entry("isotropic hardening parameter","0.", + Patterns::Double(), + "Isotropic hardening parameter of the material"); + prm.declare_entry("show stresses", "false", + Patterns::Bool(), + "Whether illustrates the stresses and von Mises stresses or not."); + + + } + + + // @sect4{The PlasticityContactProblem constructor} + + // Given the declarations of member variables as well as the + // declarations of run-time parameters that are read from the input + // file, there is nothing surprising in this constructor. In the body + // we initialize the mesh refinement strategy and the output directory, + // creating such a directory if necessary. + template + ElastoPlasticProblem:: + ElastoPlasticProblem (const ParameterHandler &prm) + : + mpi_communicator(MPI_COMM_WORLD), + n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)), + this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)), + pcout(std::cout, this_mpi_process == 0), + computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never, + TimerOutput::wall_times), + + n_initial_global_refinements (prm.get_integer("number of initial refinements")), + triangulation(mpi_communicator), + fe_degree (prm.get_integer("polynomial degree")), + fe(FE_Q(QGaussLobatto<1>(fe_degree+1)), dim), + dof_handler(triangulation), + quadrature_formula (fe_degree + 1), + face_quadrature_formula (fe_degree + 1), + + e_modulus (prm.get_double("elasticity modulus")), + nu (prm.get_double("Poissons ratio")), + sigma_0(prm.get_double("yield stress")), + gamma (prm.get_double("isotropic hardening parameter")), + constitutive_law (e_modulus, + nu, + sigma_0, + gamma), + + base_mesh (prm.get("base mesh")), + + transfer_solution (prm.get_bool("transfer solution")), + table_results(), + table_results_2(), + table_results_3(), + max_relative_error (prm.get_double("maximum relative error")), + show_stresses (prm.get_bool("show stresses")) + { + std::string strat = prm.get("refinement strategy"); + if (strat == "global") + refinement_strategy = RefinementStrategy::refine_global; + else if (strat == "percentage") + refinement_strategy = RefinementStrategy::refine_percentage; + else + AssertThrow (false, ExcNotImplemented()); + + strat = prm.get("error estimation strategy"); + if (strat == "kelly_error") + error_estimation_strategy = ErrorEstimationStrategy::kelly_error; + else if (strat == "residual_error") + error_estimation_strategy = ErrorEstimationStrategy::residual_error; + else if (strat == "weighted_residual_error") + error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error; + else + AssertThrow(false, ExcNotImplemented()); + + output_dir = prm.get("output directory"); + if (output_dir != "" && *(output_dir.rbegin()) != '/') + output_dir += "/"; + mkdir(output_dir.c_str(), 0777); + + pcout << " Using output directory '" << output_dir << "'" << std::endl; + pcout << " FE degree " << fe_degree << std::endl; + pcout << " transfer solution " + << (transfer_solution ? "true" : "false") << std::endl; + } + + + + // @sect4{PlasticityContactProblem::make_grid} + + // The next block deals with constructing the starting mesh. + // We will use the following helper function and the first + // block of the make_grid() to construct a + // mesh that corresponds to a half sphere. deal.II has a function + // that creates such a mesh, but it is in the wrong location + // and facing the wrong direction, so we need to shift and rotate + // it a bit before using it. + // + // For later reference, as described in the documentation of + // GridGenerator::half_hyper_ball(), the flat surface of the halfsphere + // has boundary indicator zero, while the remainder has boundary + // indicator one. + Point<3> + rotate_half_sphere (const Point<3> &in) + { + return Point<3>(in(2), in(1), -in(0)); + } + + template + void + ElastoPlasticProblem::make_grid () + { + if (base_mesh == "Timoshenko beam") + { + AssertThrow (dim == 2, ExcNotImplemented()); + + const double length = .48, + depth = .12; + + const Point point_1(0, -depth/2), + point_2(length, depth/2); + + std::vector repetitions(2); + repetitions[0] = 4; + repetitions[1] = 1; + GridGenerator::subdivided_hyper_rectangle(triangulation, repetitions, point_1, point_2); + + + // give the indicators to boundaries for specification, + // + // ________100______ + // | | + // 0 | | 5 + // |________________| + // 100 + // 0 to essential boundary conditions (left edge) which are as default + // 100 to the null boundaries (upper and lower edges) where we do not need to take care of them + // 5 to the natural boundaries (right edge) for imposing the traction force + typename Triangulation::cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + { + for (unsigned int face=0; face!=GeometryInfo::faces_per_cell; ++face) + { + if ( std::fabs(cell->face(face)->center()(0)-length) < 1e-12 ) + { + cell->face(face)->set_boundary_indicator(5); + }else if ( ( std::fabs(cell->face(face)->center()(1)-(depth/2)) < 1e-12 ) + || + ( std::fabs(cell->face(face)->center()(1)-(-depth/2)) < 1e-12 ) ) + { + cell->face(face)->set_boundary_indicator(100); + } + + } + } + + triangulation.refine_global(n_initial_global_refinements); + + }else if (base_mesh == "Thick_tube_internal_pressure") + { + // Example 1 from the paper: Zhong Z., .... A new numerical method for determining + // collapse load-carrying capacity of structure made of elasto-plastic material, + // J. Cent. South Univ. (2014) 21: 398-404 + AssertThrow (dim == 2, ExcNotImplemented()); + + const Point center(0, 0); + const double inner_radius = .1, + outer_radius = .2; + GridGenerator::quarter_hyper_shell(triangulation, + center, inner_radius, outer_radius, + 0, true); + + // give the indicators to boundaries for specification, + + /* _____ + | \ + | \ + 2 | \ 1 + |_ \ + \ \ + 0 \ | + |________| + 3 + */ + // 0 - inner boundary - natural boundary condition - impose the traction force + // 1 - outer boundary - free boundary - we do not need to take care of them + // 2 - left boundary - essential boundary condition - constrained to move along the x direction + // 3 - bottom boundary - essential boundary condition - constrained to move along the y direction + + const HyperBallBoundary inner_boundary_description(center, inner_radius); + triangulation.set_boundary (0, inner_boundary_description); + + const HyperBallBoundary outer_boundary_description(center, outer_radius); + triangulation.set_boundary (1, outer_boundary_description); + + triangulation.refine_global(n_initial_global_refinements); + + triangulation.set_boundary (0); + triangulation.set_boundary (1); + + }else if (base_mesh == "Perforated_strip_tension") + { + // Example 2 from the paper: Zhong Z., .... A new numerical method for determining + // collapse load-carrying capacity of structure made of elasto-plastic material, + // J. Cent. South Univ. (2014) 21: 398-404 + AssertThrow (dim == 3, ExcNotImplemented()); + + const int dim_2d = 2; + const Point center_2d(0, 0); + const double inner_radius = 0.05, + outer_radius = 0.1, + height = 0.18, + thickness = 0.004; +// thickness = 0.01; + + Triangulation triangulation_1, + triangulation_2, + triangulation_2d; + + const double eps = 1e-7 * inner_radius; + { + Point point; + + GridGenerator::quarter_hyper_shell(triangulation_1, + center_2d, inner_radius, outer_radius, + 2); + + // Modify the triangulation_1 + typename Triangulation::active_cell_iterator + cell = triangulation_1.begin_active(), + endc = triangulation_1.end(); + std::vector treated_vertices(triangulation_1.n_vertices(), false); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary() && cell->face(f)->center()(0)>eps && + cell->face(f)->center()(1)>eps ) + { + // distance of the face center from the center + point(0) = cell->face(f)->center()(0) - center_2d(0); + point(1) = cell->face(f)->center()(1) - center_2d(1); + if ( point.norm() > (inner_radius + eps) ) + { + for (unsigned int v=0; v < GeometryInfo::vertices_per_face; ++v) + { + unsigned int vv = cell->face(f)->vertex_index(v); + if (treated_vertices[vv] == false) + { + treated_vertices[vv] = true; + if (vv==1) + { + cell->face(f)->vertex(v) = center_2d+Point(outer_radius,outer_radius); + } + } + } + } + + } + } + + } + + // Make the triangulation_2, a rectangular above the triangulation_1 + { + const Point point1 (0, outer_radius), + point2 (outer_radius, height); + + GridGenerator::hyper_rectangle(triangulation_2, point1, point2); + + } + + // make the triangulation_2d and refine it + { + // Merge the two triangulation_1 and triangulation_2 + GridGenerator::merge_triangulations(triangulation_1, triangulation_2, triangulation_2d); + + // Assign boundary indicators to the boundary faces + /* + * + * /\ y + * | + * _____3_____ + * | | + * | | + * 4 | | + * | | + * | | 2 + * |_ | + * \ | + * 10 \ | + * |______| ____________\ x + * 1 / + */ + { + typename Triangulation::active_cell_iterator + cell = triangulation_2d.begin_active(), + endc = triangulation_2d.end(); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + if ( std::fabs(cell->face(f)->center()(1)) < eps ) + { + cell->face(f)->set_boundary_indicator(1); + }else if ( std::fabs(cell->face(f)->center()(0)-outer_radius) < eps ) + { + cell->face(f)->set_boundary_indicator(2); + }else if ( std::fabs(cell->face(f)->center()(1)-height) < eps ) + { + cell->face(f)->set_boundary_indicator(3); + }else if ( std::fabs(cell->face(f)->center()(0)) < eps ) + { + cell->face(f)->set_boundary_indicator(4); + }else + { + cell->face(f)->set_all_boundary_indicators(10); + } + + } + } + } + + } + + const HyperBallBoundary inner_boundary_description(center_2d, inner_radius); + triangulation_2d.set_boundary (10, inner_boundary_description); + + triangulation_2d.refine_global(3); + + triangulation_2d.set_boundary (10); + } + + // Extrude the triangulation_2d and make it 3d +// GridGenerator::extrude_triangulation(triangulation_2d, +// 2, thickness, triangulation); + extrude_triangulation(triangulation_2d, + 2, thickness, triangulation); + + // Assign boundary indicators to the boundary faces + /* + * + * /\ y + * | + * _____3_____ + * | | + * | | + * 4 | | + * | 5|6 | + * | | 2 + * |_ | + * \ | + * 10 \ | + * |______| ____________\ x + * 1 / + */ + { + Point dist_vector; + Point center(center_2d(0), center_2d(1), 0); + + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + dist_vector = cell->face(f)->center() - center; + + if ( std::fabs(dist_vector(1)) < eps ) + { + cell->face(f)->set_boundary_indicator(1); + }else if ( std::fabs(dist_vector(0)-outer_radius) < eps ) + { + cell->face(f)->set_boundary_indicator(2); + }else if ( std::fabs(dist_vector(1)-height) < eps ) + { + cell->face(f)->set_boundary_indicator(3); + }else if ( std::fabs(dist_vector(0)) < eps ) + { + cell->face(f)->set_boundary_indicator(4); + }else if ( std::fabs(dist_vector(2)) < eps ) + { + cell->face(f)->set_boundary_indicator(5); + }else if ( std::fabs(dist_vector(2)-thickness) < eps ) + { + cell->face(f)->set_boundary_indicator(6); + }else + { + cell->face(f)->set_all_boundary_indicators(10); + } + + } + } + } + + } + + const CylinderBoundary inner_boundary_description(inner_radius, 2); + triangulation.set_boundary (10, inner_boundary_description); + + triangulation.refine_global(n_initial_global_refinements); + + triangulation.set_boundary (10); + + }else if (base_mesh == "Cantiliver_beam_3d") + { + // A rectangular tube made of Aluminium + // http://www.google.de/imgres?imgurl=http%3A%2F%2Fwww.americanaluminum.com%2Fimages%2Fstockshape-rectangletube.gif&imgrefurl=http%3A%2F%2Fwww.americanaluminum.com%2Fstandard%2FrectangleTube&h=280&w=300&tbnid=VPDNh4-DJz4wyM%3A&zoom=1&docid=9DoGJCkOeFqiSM&ei=L1AuVfG5GMvtO7DggdAF&tbm=isch&client=ubuntu&iact=rc&uact=3&dur=419&page=1&start=0&ndsp=33&ved=0CGYQrQMwFQ + // approximation of beam 17250 + // units are in meter + + AssertThrow (dim == 3, ExcNotImplemented()); + + const int dim_2d = 2; + + const double length = .7, + width = 80e-3, + height = 200e-3, + thickness_web = 10e-3, + thickness_flange = 10e-3; + + Triangulation triangulation_b, + triangulation_t, + triangulation_l, + triangulation_r, + triangulation_2d; + + const double eps = 1e-7 * width; + // Make the triangulation_b, a rectangular at the bottom of rectangular tube + { + const Point point1 (-width/2, -height/2), + point2 (width/2, -(height/2)+thickness_flange); + + std::vector repetitions(dim_2d); + repetitions[0] = 8; + repetitions[1] = 1; + + GridGenerator::subdivided_hyper_rectangle(triangulation_b, repetitions, point1, point2); + } + + // Make the triangulation_t, a rectangular at the top of rectangular tube + { + const Point point1 (-width/2, (height/2)-thickness_flange), + point2 (width/2, height/2); + + std::vector repetitions(dim_2d); + repetitions[0] = 8; + repetitions[1] = 1; + + GridGenerator::subdivided_hyper_rectangle(triangulation_t, repetitions, point1, point2); + } + + // Make the triangulation_l, a rectangular at the left of rectangular tube + { + const Point point1 (-width/2, -(height/2)+thickness_flange), + point2 (-(width/2)+thickness_web, (height/2)-thickness_flange); + + std::vector repetitions(dim_2d); + repetitions[0] = 1; + repetitions[1] = 18; + + GridGenerator::subdivided_hyper_rectangle(triangulation_l, repetitions, point1, point2); + } + + // Make the triangulation_r, a rectangular at the right of rectangular tube + { + const Point point1 ((width/2)-thickness_web, -(height/2)+thickness_flange), + point2 (width/2, (height/2)-thickness_flange); + + std::vector repetitions(dim_2d); + repetitions[0] = 1; + repetitions[1] = 18; + + GridGenerator::subdivided_hyper_rectangle(triangulation_r, repetitions, point1, point2); + } + + // make the triangulation_2d + { + // merging every two triangles to make triangulation_2d + Triangulation triangulation_bl, + triangulation_blr; + + GridGenerator::merge_triangulations(triangulation_b, triangulation_l, triangulation_bl); + GridGenerator::merge_triangulations(triangulation_bl, triangulation_r, triangulation_blr); + GridGenerator::merge_triangulations(triangulation_blr, triangulation_t, triangulation_2d); + } + + // Extrude the triangulation_2d and make it 3d + const unsigned int n_slices = length*1000/20 + 1; + extrude_triangulation(triangulation_2d, + n_slices, length, triangulation); + + // Assign boundary indicators to the boundary faces + /* + * + * A + * ---------*---------- + * / /| + * / / | + * / / | + * / 2 length / | + * / / | + * / / | + * / / | + * / width / | + * -------------------- | + * | --------1-------. | | + * | : : | | + * | : : |h | + * | : y z : |e | + * | : | / : |i / + * |1: |___ x :1|g / + * | : : |h / + * | : : |t / + * | : : | / + * | : : | / + * | ----------------- |/ + * ---------1----------/ + * + * face id: + * Essential boundary condition: + * 1: z = 0: clamped, fixed in x, y and z directions + * Natural/Newmann boundary condition: + * 2: y = height/2: traction face: pressure on the surface + * Quantity of interest: + * displacement at Point A (x=0, y=height/2, z=length) + */ + { + Point dist_vector; + Point center(0, 0, 0); + + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + dist_vector = cell->face(f)->center() - center; + + if ( std::fabs(dist_vector(2)) < eps ) + { + cell->face(f)->set_boundary_indicator(1); + }else if ( std::fabs(dist_vector(1)-(height/2)) < eps ) + { + cell->face(f)->set_boundary_indicator(2); + }else + { + cell->face(f)->set_all_boundary_indicators(0); + } + + } + } + } + + } + + triangulation.refine_global(n_initial_global_refinements); + + }else + { + AssertThrow(false, ExcNotImplemented()); + } + + pcout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + } + + + + // @sect4{PlasticityContactProblem::setup_system} + + // The next piece in the puzzle is to set up the DoFHandler, resize + // vectors and take care of various other status variables such as + // index sets and constraint matrices. + // + // In the following, each group of operations is put into a brace-enclosed + // block that is being timed by the variable declared at the top of the + // block (the constructor of the TimerOutput::Scope variable starts the + // timed section, the destructor that is called at the end of the block + // stops it again). + template + void + ElastoPlasticProblem::setup_system () + { + /* setup dofs and get index sets for locally owned and relevant dofs */ + TimerOutput::Scope t(computing_timer, "Setup"); + { + TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs"); + dof_handler.distribute_dofs(fe); + pcout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs(dof_handler, + locally_relevant_dofs); + } + + /* setup hanging nodes and Dirichlet constraints */ + { + TimerOutput::Scope t(computing_timer, "Setup: constraints"); + constraints_hanging_nodes.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, + constraints_hanging_nodes); + constraints_hanging_nodes.close(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + compute_dirichlet_constraints(); + } + + /* initialization of vectors*/ + { + TimerOutput::Scope t(computing_timer, "Setup: vectors"); + if (timestep_no==1 || current_refinement_cycle!=0) + { + solution.reinit(locally_relevant_dofs, mpi_communicator); + } + incremental_displacement.reinit(locally_relevant_dofs, mpi_communicator); + newton_rhs.reinit(locally_owned_dofs, mpi_communicator); + newton_rhs_residual.reinit(locally_owned_dofs, mpi_communicator); + fraction_of_plastic_q_points_per_cell.reinit(triangulation.n_active_cells()); + } + + // Finally, we set up sparsity patterns and matrices. + // We temporarily (ab)use the system matrix to also build the (diagonal) + // matrix that we use in eliminating degrees of freedom that are in contact + // with the obstacle, but we then immediately set the Newton matrix back + // to zero. + { + TimerOutput::Scope t(computing_timer, "Setup: matrix"); + TrilinosWrappers::SparsityPattern sp(locally_owned_dofs, + mpi_communicator); + + DoFTools::make_sparsity_pattern(dof_handler, sp, + constraints_dirichlet_and_hanging_nodes, false, + this_mpi_process); + sp.compress(); + newton_matrix.reinit(sp); + } + } + + + // @sect4{PlasticityContactProblem::compute_dirichlet_constraints} + + // This function, broken out of the preceding one, computes the constraints + // associated with Dirichlet-type boundary conditions and puts them into the + // constraints_dirichlet_and_hanging_nodes variable by merging + // with the constraints that come from hanging nodes. + // + // As laid out in the introduction, we need to distinguish between two + // cases: + // - If the domain is a box, we set the displacement to zero at the bottom, + // and allow vertical movement in z-direction along the sides. As + // shown in the make_grid() function, the former corresponds + // to boundary indicator 6, the latter to 8. + // - If the domain is a half sphere, then we impose zero displacement along + // the curved part of the boundary, associated with boundary indicator zero. + template + void + ElastoPlasticProblem::compute_dirichlet_constraints () + { + constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs); + constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes); + + std::vector component_mask(dim); + + if (base_mesh == "Timoshenko beam") + { + VectorTools::interpolate_boundary_values(dof_handler, + 0, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes, + ComponentMask()); + }else if (base_mesh == "Thick_tube_internal_pressure") + { + // the boundary x = 0 + component_mask[0] = true; component_mask[1] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 2, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes, + component_mask); + // the boundary y = 0 + component_mask[0] = false; component_mask[1] = true; + VectorTools::interpolate_boundary_values (dof_handler, + 3, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes, + component_mask); + }else if (base_mesh == "Perforated_strip_tension") + { + // the boundary x = 0 + component_mask[0] = true; component_mask[1] = false; component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 4, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes, + component_mask); + // the boundary y = 0 + component_mask[0] = false; component_mask[1] = true; component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 1, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes, + component_mask); + // the boundary y = imposed incremental displacement + component_mask[0] = false; component_mask[1] = true; component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 3, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes, + component_mask); + }else if (base_mesh == "Cantiliver_beam_3d") + { + // the boundary x = y = z = 0 + component_mask[0] = true; component_mask[1] = true; component_mask[2] = true; + VectorTools::interpolate_boundary_values (dof_handler, + 1, + EquationData::IncrementalBoundaryValues(present_time, end_time), + constraints_dirichlet_and_hanging_nodes, + component_mask); + }else + { + AssertThrow(false, ExcNotImplemented()); + } + + + constraints_dirichlet_and_hanging_nodes.close(); + } + + + // @sect4{PlasticityContactProblem::assemble_newton_system} + + // Given the complexity of the problem, it may come as a bit of a surprise + // that assembling the linear system we have to solve in each Newton iteration + // is actually fairly straightforward. The following function builds the Newton + // right hand side and Newton matrix. It looks fairly innocent because the + // heavy lifting happens in the call to + // ConstitutiveLaw::get_linearized_stress_strain_tensors() and in + // particular in ConstraintMatrix::distribute_local_to_global(), using the + // constraints we have previously computed. + template + void + ElastoPlasticProblem:: + assemble_newton_system (const TrilinosWrappers::MPI::Vector &linearization_point, + const TrilinosWrappers::MPI::Vector &delta_linearization_point) + { + TimerOutput::Scope t(computing_timer, "Assembling"); + + types::boundary_id traction_surface_id; + if (base_mesh == "Timoshenko beam") + { + traction_surface_id = 5; + }else if (base_mesh == "Thick_tube_internal_pressure") + { + traction_surface_id = 0; + }else if (base_mesh == "Cantiliver_beam_3d") + { + traction_surface_id = 2; + } + + FEValues fe_values(fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + + const EquationData::BodyForce body_force; + std::vector > body_force_values(n_q_points, + Vector(dim)); + + const EquationData:: + IncrementalBoundaryForce boundary_force(present_time, end_time); + std::vector > boundary_force_values(n_face_q_points, + Vector(dim)); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + +// std::vector > strain_tensor(n_q_points); + std::vector > incremental_strain_tensor(n_q_points); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell != endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + + fe_values[displacement].get_function_symmetric_gradients(delta_linearization_point, + incremental_strain_tensor); + + // For assembling the local right hand side contributions, we need + // to access the prior linearized stress value in this quadrature + // point. To get it, we use the user pointer of this cell that + // points into the global array to the quadrature point data + // corresponding to the first quadrature point of the present cell, + // and then add an offset corresponding to the index of the + // quadrature point we presently consider: + const PointHistory *local_quadrature_points_history + = reinterpret_cast*>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + + // In addition, we need the values of the external body forces at + // the quadrature points on this cell: + body_force.vector_value_list(fe_values.get_quadrature_points(), + body_force_values); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + SymmetricTensor<2, dim> tmp_strain_tensor_qpoint; + tmp_strain_tensor_qpoint = local_quadrature_points_history[q_point].old_strain + + incremental_strain_tensor[q_point]; + + SymmetricTensor<4, dim> stress_strain_tensor_linearized; + SymmetricTensor<4, dim> stress_strain_tensor; + constitutive_law.get_linearized_stress_strain_tensors(tmp_strain_tensor_qpoint, + stress_strain_tensor_linearized, + stress_strain_tensor); + + Tensor<1, dim> rhs_values_body_force; + for (unsigned int i = 0; i < dim; ++i) + { + rhs_values_body_force[i] = body_force_values[q_point][i]; + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + // Having computed the stress-strain tensor and its linearization, + // we can now put together the parts of the matrix and right hand side. + // In both, we need the linearized stress-strain tensor times the + // symmetric gradient of $\varphi_i$, i.e. the term $I_\Pi\varepsilon(\varphi_i)$, + // so we introduce an abbreviation of this term. Recall that the + // matrix corresponds to the bilinear form + // $A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))$ in the + // notation of the accompanying publication, whereas the right + // hand side is $F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))$ + // where $u$ is the current linearization points (typically the last solution). + // This might suggest that the right hand side will be zero if the material + // is completely elastic (where $I_\Pi=P_\Pi$) but this ignores the fact + // that the right hand side will also contain contributions from + // non-homogeneous constraints due to the contact. + // + // The code block that follows this adds contributions that are due to + // boundary forces, should there be any. + const SymmetricTensor<2, dim> + stress_phi_i = stress_strain_tensor_linearized + * fe_values[displacement].symmetric_gradient(i, q_point); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += (stress_phi_i + * fe_values[displacement].symmetric_gradient(j, q_point) + * fe_values.JxW(q_point)); + + cell_rhs(i) += ( + ( stress_phi_i + * incremental_strain_tensor[q_point] ) + - + ( ( stress_strain_tensor + * fe_values[displacement].symmetric_gradient(i, q_point)) + * tmp_strain_tensor_qpoint ) + + + ( fe_values[displacement].value(i, q_point) + * rhs_values_body_force ) + ) * fe_values.JxW(q_point); + + } + } + + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary() + && + cell->face(face)->boundary_indicator() == traction_surface_id) + { + fe_values_face.reinit(cell, face); + + boundary_force.vector_value_list(fe_values_face.get_quadrature_points(), + boundary_force_values); + + for (unsigned int q_point=0; q_point rhs_values; + for (unsigned int i = 0; i < dim; ++i) + { + rhs_values[i] = boundary_force_values[q_point][i]; + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) + * rhs_values + * fe_values_face.JxW(q_point)); + } + } + + cell->get_dof_indices(local_dof_indices); + constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_matrix, cell_rhs, + local_dof_indices, + newton_matrix, + newton_rhs, + true); + + } + + newton_matrix.compress(VectorOperation::add); + newton_rhs.compress(VectorOperation::add); + } + + + + // @sect4{PlasticityContactProblem::compute_nonlinear_residual} + + // The following function computes the nonlinear residual of the equation + // given the current solution (or any other linearization point). This + // is needed in the linear search algorithm where we need to try various + // linear combinations of previous and current (trial) solution to + // compute the (real, globalized) solution of the current Newton step. + // + // That said, in a slight abuse of the name of the function, it actually + // does significantly more. For example, it also computes the vector + // that corresponds to the Newton residual but without eliminating + // constrained degrees of freedom. We need this vector to compute contact + // forces and, ultimately, to compute the next active set. Likewise, by + // keeping track of how many quadrature points we encounter on each cell + // that show plastic yielding, we also compute the + // fraction_of_plastic_q_points_per_cell vector that we + // can later output to visualize the plastic zone. In both of these cases, + // the results are not necessary as part of the line search, and so we may + // be wasting a small amount of time computing them. At the same time, this + // information appears as a natural by-product of what we need to do here + // anyway, and we want to collect it once at the end of each Newton + // step, so we may as well do it here. + // + // The actual implementation of this function should be rather obvious: + template + void + ElastoPlasticProblem:: + compute_nonlinear_residual (const TrilinosWrappers::MPI::Vector &linearization_point) + { + types::boundary_id traction_surface_id; + if (base_mesh == "Timoshenko beam") + { + traction_surface_id = 5; + }else if (base_mesh == "Thick_tube_internal_pressure") + { + traction_surface_id = 0; + }else if (base_mesh == "Cantiliver_beam_3d") + { + traction_surface_id = 2; + } + + FEValues fe_values(fe, quadrature_formula, + update_values | update_gradients | update_quadrature_points | + update_JxW_values); + + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const EquationData::BodyForce body_force; + std::vector > body_force_values(n_q_points, + Vector(dim)); + + const EquationData:: + IncrementalBoundaryForce boundary_force(present_time, end_time); + std::vector > boundary_force_values(n_face_q_points, + Vector(dim)); + + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const FEValuesExtractors::Vector displacement(0); + + newton_rhs_residual = 0; + + fraction_of_plastic_q_points_per_cell = 0; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + unsigned int cell_number = 0; + for (; cell != endc; ++cell, ++cell_number) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_rhs = 0; + + std::vector > strain_tensors(n_q_points); + fe_values[displacement].get_function_symmetric_gradients(linearization_point, + strain_tensors); + + body_force.vector_value_list(fe_values.get_quadrature_points(), + body_force_values); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + SymmetricTensor<4, dim> stress_strain_tensor; + const bool q_point_is_plastic + = constitutive_law.get_stress_strain_tensor(strain_tensors[q_point], + stress_strain_tensor); + if (q_point_is_plastic) + ++fraction_of_plastic_q_points_per_cell(cell_number); + + Tensor<1, dim> rhs_values_body_force; + for (unsigned int i = 0; i < dim; ++i) + { + rhs_values_body_force[i] = body_force_values[q_point][i]; + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + cell_rhs(i) += (fe_values[displacement].value(i, q_point) + * rhs_values_body_force + - + strain_tensors[q_point] + * stress_strain_tensor + * fe_values[displacement].symmetric_gradient(i, q_point) + ) + * fe_values.JxW(q_point); + + Tensor<1, dim> rhs_values; + rhs_values = 0; + cell_rhs(i) += (fe_values[displacement].value(i, q_point) + * rhs_values + * fe_values.JxW(q_point)); + } + } + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face(face)->at_boundary() + && cell->face(face)->boundary_indicator() == traction_surface_id) + { + fe_values_face.reinit(cell, face); + + boundary_force.vector_value_list(fe_values_face.get_quadrature_points(), + boundary_force_values); + + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) + { + Tensor<1, dim> rhs_values; + for (unsigned int i = 0; i < dim; ++i) + { + rhs_values[i] = boundary_force_values[q_point][i]; + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values + * fe_values_face.JxW(q_point)); + } + } + + cell->get_dof_indices(local_dof_indices); + constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_rhs, + local_dof_indices, + newton_rhs_residual); + + } + + fraction_of_plastic_q_points_per_cell /= quadrature_formula.size(); + newton_rhs_residual.compress(VectorOperation::add); + + } + + + + + + // @sect4{PlasticityContactProblem::solve_newton_system} + + // The last piece before we can discuss the actual Newton iteration + // on a single mesh is the solver for the linear systems. There are + // a couple of complications that slightly obscure the code, but + // mostly it is just setup then solve. Among the complications are: + // + // - For the hanging nodes we have to apply + // the ConstraintMatrix::set_zero function to newton_rhs. + // This is necessary if a hanging node with solution value $x_0$ + // has one neighbor with value $x_1$ which is in contact with the + // obstacle and one neighbor $x_2$ which is not in contact. Because + // the update for the former will be prescribed, the hanging node constraint + // will have an inhomogeneity and will look like $x_0 = x_1/2 + \text{gap}/2$. + // So the corresponding entries in the + // ride-hang-side are non-zero with a + // meaningless value. These values we have to + // to set to zero. + // - Like in step-40, we need to shuffle between vectors that do and do + // do not have ghost elements when solving or using the solution. + // + // The rest of the function is similar to step-40 and + // step-41 except that we use a BiCGStab solver + // instead of CG. This is due to the fact that for very small hardening + // parameters $\gamma$, the linear system becomes almost semidefinite though + // still symmetric. BiCGStab appears to have an easier time with such linear + // systems. + template + void + ElastoPlasticProblem::solve_newton_system () + { + TimerOutput::Scope t(computing_timer, "Solve"); + + TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs, mpi_communicator); + distributed_solution = incremental_displacement; + + constraints_hanging_nodes.set_zero(distributed_solution); + constraints_hanging_nodes.set_zero(newton_rhs); + + // ------- Solver Bicgstab --- Preconditioner AMG ------------------- +// TrilinosWrappers::PreconditionAMG preconditioner; +// { +// TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner"); +// +// std::vector > constant_modes; +// DoFTools::extract_constant_modes(dof_handler, ComponentMask(), +// constant_modes); +// +// TrilinosWrappers::PreconditionAMG::AdditionalData additional_data; +// additional_data.constant_modes = constant_modes; +// additional_data.elliptic = true; +// additional_data.n_cycles = 1; +// additional_data.w_cycle = false; +// additional_data.output_details = false; +// additional_data.smoother_sweeps = 2; +// additional_data.aggregation_threshold = 1e-2; +// +// preconditioner.initialize(newton_matrix, additional_data); +// } + +// { +// TimerOutput::Scope t(computing_timer, "Solve: iterate"); +// +// TrilinosWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator); +// +//// const double relative_accuracy = 1e-8; +// const double relative_accuracy = 1e-2; +// const double solver_tolerance = relative_accuracy +// * newton_matrix.residual(tmp, distributed_solution, +// newton_rhs); +// +// SolverControl solver_control(newton_matrix.m(), +// solver_tolerance); +// SolverBicgstab solver(solver_control); +// solver.solve(newton_matrix, distributed_solution, +// newton_rhs, preconditioner); +// +// pcout << " Error: " << solver_control.initial_value() +// << " -> " << solver_control.last_value() << " in " +// << solver_control.last_step() << " Bicgstab iterations." +// << std::endl; +// } + + // ------- Solver CG --- Preconditioner SSOR ------------------- + TrilinosWrappers::PreconditionSSOR preconditioner; + { + TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner"); + + TrilinosWrappers::PreconditionSSOR::AdditionalData additional_data; + preconditioner.initialize(newton_matrix, additional_data); + } + + { + TimerOutput::Scope t(computing_timer, "Solve: iterate"); + + TrilinosWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator); + +// const double relative_accuracy = 1e-8; + const double relative_accuracy = 1e-2; + const double solver_tolerance = relative_accuracy + * newton_matrix.residual(tmp, distributed_solution, + newton_rhs); + +// SolverControl solver_control(newton_matrix.m(), +// solver_tolerance); + SolverControl solver_control(10*newton_matrix.m(), + solver_tolerance); + SolverCG solver(solver_control); + solver.solve(newton_matrix, distributed_solution, + newton_rhs, preconditioner); + + pcout << " Error: " << solver_control.initial_value() + << " -> " << solver_control.last_value() << " in " + << solver_control.last_step() << " CG iterations." + << std::endl; + } + // ........................................................ + + constraints_dirichlet_and_hanging_nodes.distribute(distributed_solution); + + incremental_displacement = distributed_solution; + } + + + // @sect4{PlasticityContactProblem::solve_newton} + + // This is, finally, the function that implements the damped Newton method + // on the current mesh. There are two nested loops: the outer loop for the Newton + // iteration and the inner loop for the line search which + // will be used only if necessary. To obtain a good and reasonable + // starting value we solve an elastic problem in very first Newton step on each + // mesh (or only on the first mesh if we transfer solutions between meshes). We + // do so by setting the yield stress to an unreasonably large value in these + // iterations and then setting it back to the correct value in subsequent + // iterations. + // + // Other than this, the top part of this function should be reasonably + // obvious: + template + void + ElastoPlasticProblem::solve_newton () + { + TrilinosWrappers::MPI::Vector old_solution(locally_owned_dofs, mpi_communicator); + TrilinosWrappers::MPI::Vector residual(locally_owned_dofs, mpi_communicator); + TrilinosWrappers::MPI::Vector tmp_vector(locally_owned_dofs, mpi_communicator); + TrilinosWrappers::MPI::Vector locally_relevant_tmp_vector(locally_relevant_dofs, mpi_communicator); + TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs, mpi_communicator); + TrilinosWrappers::MPI::Vector tmp_solution(locally_owned_dofs, mpi_communicator); + + double residual_norm; + double previous_residual_norm = -std::numeric_limits::max(); + + double disp_norm, + previous_disp_norm = 0; + + const double correct_sigma = sigma_0; + + const unsigned int max_newton_iter = 100; + + for (unsigned int newton_step = 1; newton_step <= max_newton_iter; ++newton_step) + { + if (newton_step == 1 + && + ((transfer_solution && timestep_no == 1) + || + !transfer_solution)) + constitutive_law.set_sigma_0(1e+10); + else + constitutive_law.set_sigma_0(correct_sigma); + + pcout << " " << std::endl; + pcout << " Newton iteration " << newton_step << std::endl; + + pcout << " Assembling system... " << std::endl; + newton_matrix = 0; + newton_rhs = 0; + newton_rhs_residual = 0; + + tmp_solution = solution; + tmp_solution += incremental_displacement; + assemble_newton_system(tmp_solution, + incremental_displacement); + + pcout << " Solving system... " << std::endl; + solve_newton_system(); + + // It gets a bit more hairy after we have computed the + // trial solution $\tilde{\mathbf u}$ of the current Newton step. + // We handle a highly nonlinear problem so we have to damp + // Newton's method using a line search. To understand how we do this, + // recall that in our formulation, we compute a trial solution + // in each Newton step and not the update between old and new solution. + // Since the solution set is a convex set, we will use a line + // search that tries linear combinations of the + // previous and the trial solution to guarantee that the + // damped solution is in our solution set again. + // At most we apply 5 damping steps. + // + // There are exceptions to when we use a line search. First, + // if this is the first Newton step on any mesh, then we don't have + // any point to compare the residual to, so we always accept a full + // step. Likewise, if this is the second Newton step on the first mesh (or + // the second on any mesh if we don't transfer solutions from + // mesh to mesh), then we have computed the first of these steps using + // just an elastic model (see how we set the yield stress sigma to + // an unreasonably large value above). In this case, the first Newton + // solution was a purely elastic one, the second one a plastic one, + // and any linear combination would not necessarily be expected to + // lie in the feasible set -- so we just accept the solution we just + // got. + // + // In either of these two cases, we bypass the line search and just + // update residual and other vectors as necessary. + if ((newton_step==1) + || + (transfer_solution && newton_step == 2 && current_refinement_cycle == 0) + || + (!transfer_solution && newton_step == 2)) + { + tmp_solution = solution; + tmp_solution += incremental_displacement; + compute_nonlinear_residual(tmp_solution); + old_solution = incremental_displacement; + + residual = newton_rhs_residual; + + residual.compress(VectorOperation::insert); + + residual_norm = residual.l2_norm(); + + pcout << " Accepting Newton solution with residual: " + << residual_norm << std::endl; + } + else + { + for (unsigned int i = 0; i < 5; i++) + { + distributed_solution = incremental_displacement; + + const double alpha = std::pow(0.5, static_cast(i)); + tmp_vector = old_solution; + tmp_vector.sadd(1 - alpha, alpha, distributed_solution); + + TimerOutput::Scope t(computing_timer, "Residual and lambda"); + + locally_relevant_tmp_vector = tmp_vector; + tmp_solution = solution; + tmp_solution += locally_relevant_tmp_vector; + compute_nonlinear_residual(tmp_solution); + residual = newton_rhs_residual; + + residual.compress(VectorOperation::insert); + + residual_norm = residual.l2_norm(); + + pcout << " Residual of the system: " + << residual_norm << std::endl + << " with a damping parameter alpha = " << alpha + << std::endl; + + if (residual_norm < previous_residual_norm) + break; + } + + incremental_displacement = tmp_vector; + old_solution = incremental_displacement; + } + + disp_norm = incremental_displacement.l2_norm(); + + + // The final step is to check for convergence. If the residual is + // less than a threshold of $10^{-10}$, then we terminate + // the iteration on the current mesh: +// if (residual_norm < 1e-10) + if (residual_norm < 1e-7) + break; + + pcout << " difference of two consecutive incremental displacement l2 norm : " + << std::abs(disp_norm - previous_disp_norm) << std::endl; + if ( std::abs(disp_norm - previous_disp_norm) < 1e-10 && + (residual_norm < 1e-5 || std::abs(residual_norm - previous_residual_norm)<1e-9) ) + { + pcout << " Convergence by difference of two consecutive solution! " << std::endl; + break; + } + + + previous_residual_norm = residual_norm; + previous_disp_norm = disp_norm; + } + } + + // @sect4{PlasticityContactProblem::compute_error} + + template + void + ElastoPlasticProblem::compute_error () + { + TrilinosWrappers::MPI::Vector tmp_solution(locally_owned_dofs, mpi_communicator); + tmp_solution = solution; + tmp_solution += incremental_displacement; + + estimated_error_per_cell.reinit (triangulation.n_active_cells()); + if (error_estimation_strategy == ErrorEstimationStrategy::kelly_error) + { + KellyErrorEstimator::estimate(dof_handler, + QGauss(fe.degree + 2), + typename FunctionMap::type(), + tmp_solution, + estimated_error_per_cell); + + }else if (error_estimation_strategy == ErrorEstimationStrategy::residual_error) + { + compute_error_residual(tmp_solution); + + }else if (error_estimation_strategy == ErrorEstimationStrategy::weighted_residual_error) + { + // make a non-parallel copy of tmp_solution + Vector copy_solution(tmp_solution); + + // the dual function definition (it should be defined previously, e.g. input file) + if (base_mesh == "Timoshenko beam") + { + double length = .48, + depth = .12; + + const Point evaluation_point(length, -depth/2); + + DualFunctional::PointValuesEvaluation dual_functional(evaluation_point); + + DualSolver dual_solver(triangulation, fe, + copy_solution, + constitutive_law, dual_functional, + timestep_no, output_dir, base_mesh, + present_time, end_time); + + dual_solver.compute_error_DWR (estimated_error_per_cell); + + }else if (base_mesh == "Thick_tube_internal_pressure") + { + const unsigned int face_id = 0; + std::vector > comp_stress(dim); + for (unsigned int i=0; i!=dim; ++i) + { + comp_stress[i].resize(dim); + for (unsigned int j=0; j!=dim; ++j) + { + comp_stress[i][j] = 1; + } + } + + DualFunctional::MeanStressFace dual_functional(face_id, comp_stress); + + DualSolver dual_solver(triangulation, fe, + copy_solution, + constitutive_law, dual_functional, + timestep_no, output_dir, base_mesh, + present_time, end_time); + + dual_solver.compute_error_DWR (estimated_error_per_cell); + + }else if (base_mesh == "Perforated_strip_tension") + { + // ......................................... + // Mean stress_yy over the bottom boundary + const unsigned int face_id = 1; + std::vector > comp_stress(dim); + for (unsigned int i=0; i!=dim; ++i) + { + comp_stress[i].resize(dim); + for (unsigned int j=0; j!=dim; ++j) + { + comp_stress[i][j] = 0; + } + } + comp_stress[1][1] = 1; + + DualFunctional::MeanStressFace dual_functional(face_id, comp_stress); + + // ......................................... + + DualSolver dual_solver(triangulation, fe, + copy_solution, + constitutive_law, dual_functional, + timestep_no, output_dir, base_mesh, + present_time, end_time); + + dual_solver.compute_error_DWR (estimated_error_per_cell); + + }else if (base_mesh == "Cantiliver_beam_3d") + { + // Quantity of interest: + // ----------------------------------------------------------- + // displacement at Point A (x=0, y=height/2, z=length) + /* + const double length = .7, + height = 200e-3; + + const Point evaluation_point(0, height/2, length); + + DualFunctional::PointValuesEvaluation dual_functional(evaluation_point); + */ + + // ----------------------------------------------------------- + // Mean stress at the specified domain is of interest. + // The interest domains are located on the bottom and top of the flanges + // close to the clamped face, z = 0 + // top domain: height/2 - thickness_flange <= y <= height/2 + // 0 <= z <= 2 * thickness_flange + // bottom domain: -height/2 <= y <= -height/2 + thickness_flange + // 0 <= z <= 2 * thickness_flange + + std::vector > comp_stress(dim); + for (unsigned int i=0; i!=dim; ++i) + { + comp_stress[i].resize(dim); + for (unsigned int j=0; j!=dim; ++j) + { + comp_stress[i][j] = 1; + } + } + DualFunctional::MeanStressDomain dual_functional(base_mesh, comp_stress); + + // ----------------------------------------------------------- + + DualSolver dual_solver(triangulation, fe, + copy_solution, + constitutive_law, dual_functional, + timestep_no, output_dir, base_mesh, + present_time, end_time); + + dual_solver.compute_error_DWR (estimated_error_per_cell); + + }else + { + AssertThrow(false, ExcNotImplemented()); + } + + + }else + { + AssertThrow(false, ExcNotImplemented()); + } + + + relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm(); + + pcout << "Estimated relative error = " << relative_error << std::endl; + + } + + template + void + ElastoPlasticProblem::compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution) + { + FEValues fe_values(fe, quadrature_formula, + update_values | + update_gradients | + update_hessians | + update_quadrature_points | + update_JxW_values); + + const unsigned int n_q_points = quadrature_formula.size(); + std::vector > strain_tensor(n_q_points); + SymmetricTensor<4, dim> stress_strain_tensor_linearized; + SymmetricTensor<4, dim> stress_strain_tensor; + Tensor<5, dim> stress_strain_tensor_grad; + std::vector > > cell_hessians (n_q_points); + for (unsigned int i=0; i!=n_q_points; ++i) + { + cell_hessians[i].resize (dim); + } + const EquationData::BodyForce body_force; + + std::vector > body_force_values (n_q_points, Vector(dim)); + const FEValuesExtractors::Vector displacement(0); + + + FEFaceValues fe_face_values_cell(fe, face_quadrature_formula, + update_values | + update_quadrature_points| + update_gradients | + update_JxW_values | + update_normal_vectors), + fe_face_values_neighbor (fe, face_quadrature_formula, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors); + FESubfaceValues fe_subface_values_cell (fe, face_quadrature_formula, + update_gradients); + + const unsigned int n_face_q_points = face_quadrature_formula.size(); + std::vector > jump_residual (n_face_q_points, Vector(dim)); + std::vector > > cell_grads(n_face_q_points); + for (unsigned int i=0; i!=n_face_q_points; ++i) + { + cell_grads[i].resize (dim); + } + std::vector > > neighbor_grads(n_face_q_points); + for (unsigned int i=0; i!=n_face_q_points; ++i) + { + neighbor_grads[i].resize (dim); + } + SymmetricTensor<2, dim> q_cell_strain_tensor; + SymmetricTensor<2, dim> q_neighbor_strain_tensor; + SymmetricTensor<4, dim> cell_stress_strain_tensor; + SymmetricTensor<4, dim> neighbor_stress_strain_tensor; + + + typename std::map::face_iterator, Vector > + face_integrals; + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + { + face_integrals[cell->face(face_no)].reinit (dim); + face_integrals[cell->face(face_no)] = -1e20; + } + } + + std::vector > error_indicators_vector; + error_indicators_vector.resize( triangulation.n_active_cells(), + Vector(dim) ); + + // ----------------- estimate_some ------------------------- + cell = dof_handler.begin_active(); + unsigned int present_cell = 0; + for (; cell!=endc; ++cell, ++present_cell) + if (cell->is_locally_owned()) + { + // --------------- integrate_over_cell ------------------- + fe_values.reinit(cell); + body_force.vector_value_list(fe_values.get_quadrature_points(), + body_force_values); + fe_values[displacement].get_function_symmetric_gradients(tmp_solution, + strain_tensor); + fe_values.get_function_hessians(tmp_solution, cell_hessians); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point], + stress_strain_tensor_linearized, + stress_strain_tensor); + constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point], + cell_hessians[q_point], + stress_strain_tensor_grad); + + for (unsigned int i=0; i!=dim; ++i) + { + error_indicators_vector[present_cell](i) += + body_force_values[q_point](i)*fe_values.JxW(q_point); + for (unsigned int j=0; j!=dim; ++j) + { + for (unsigned int k=0; k!=dim; ++k) + { + for (unsigned int l=0; l!=dim; ++l) + { + error_indicators_vector[present_cell](i) += + ( stress_strain_tensor[i][j][k][l]* + 0.5*(cell_hessians[q_point][k][l][j] + + + cell_hessians[q_point][l][k][j]) + + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l] + ) * + fe_values.JxW(q_point); + } + } + } + + } + + } + // ------------------------------------------------------- + // compute face_integrals + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + { + if (cell->face(face_no)->at_boundary()) + { + for (unsigned int id=0; id!=dim; ++id) + { + face_integrals[cell->face(face_no)](id) = 0; + } + continue; + } + + if ((cell->neighbor(face_no)->has_children() == false) && + (cell->neighbor(face_no)->level() == cell->level()) && + (cell->neighbor(face_no)->index() < cell->index())) + continue; + + if (cell->at_boundary(face_no) == false) + if (cell->neighbor(face_no)->level() < cell->level()) + continue; + + + if (cell->face(face_no)->has_children() == false) + { + // ------------- integrate_over_regular_face ----------- + fe_face_values_cell.reinit(cell, face_no); + fe_face_values_cell.get_function_grads (tmp_solution, + cell_grads); + + Assert (cell->neighbor(face_no).state() == IteratorState::valid, + ExcInternalError()); + const unsigned int + neighbor_neighbor = cell->neighbor_of_neighbor (face_no); + const typename DoFHandler::active_cell_iterator + neighbor = cell->neighbor(face_no); + + fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor); + fe_face_values_neighbor.get_function_grads (tmp_solution, + neighbor_grads); + + for (unsigned int q_point=0; q_point face_integral_vector(dim); + face_integral_vector = 0; + for (unsigned int q_point=0; q_pointface(face_no)) != face_integrals.end(), + ExcInternalError()); + + for (unsigned int i=0; i!=dim; ++i) + { + Assert (face_integrals[cell->face(face_no)](i) == -1e20, + ExcInternalError()); + face_integrals[cell->face(face_no)](i) = face_integral_vector(i); + + } + + // ----------------------------------------------------- + }else + { + // ------------- integrate_over_irregular_face --------- + const typename DoFHandler::face_iterator + face = cell->face(face_no); + const typename DoFHandler::cell_iterator + neighbor = cell->neighbor(face_no); + Assert (neighbor.state() == IteratorState::valid, + ExcInternalError()); + Assert (neighbor->has_children(), + ExcInternalError()); + + const unsigned int + neighbor_neighbor = cell->neighbor_of_neighbor (face_no); + + for (unsigned int subface_no=0; + subface_non_children(); ++subface_no) + { + const typename DoFHandler::active_cell_iterator + neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no); + Assert (neighbor_child->face(neighbor_neighbor) == + cell->face(face_no)->child(subface_no), + ExcInternalError()); + + fe_subface_values_cell.reinit (cell, face_no, subface_no); + fe_subface_values_cell.get_function_grads (tmp_solution, + cell_grads); + fe_face_values_neighbor.reinit (neighbor_child, + neighbor_neighbor); + fe_face_values_neighbor.get_function_grads (tmp_solution, + neighbor_grads); + + for (unsigned int q_point=0; q_point face_integral_vector(dim); + face_integral_vector = 0; + for (unsigned int q_point=0; q_pointface(neighbor_neighbor)](i) = face_integral_vector(i); + } + + } + + Vector sum (dim); + sum = 0; + for (unsigned int subface_no=0; + subface_non_children(); ++subface_no) + { + Assert (face_integrals.find(face->child(subface_no)) != + face_integrals.end(), + ExcInternalError()); + for (unsigned int i=0; i!=dim; ++i) + { + Assert (face_integrals[face->child(subface_no)](i) != -1e20, + ExcInternalError()); + sum(i) += face_integrals[face->child(subface_no)](i); + } + } + for (unsigned int i=0; i!=dim; ++i) + { + face_integrals[face](i) = sum(i); + } + + + // ----------------------------------------------------- + } + + + } + } + // ---------------------------------------------------------- + + present_cell=0; + cell = dof_handler.begin_active(); + for (; cell!=endc; ++cell, ++present_cell) + if (cell->is_locally_owned()) + { + for (unsigned int face_no=0; face_no::faces_per_cell; + ++face_no) + { + Assert(face_integrals.find(cell->face(face_no)) != + face_integrals.end(), + ExcInternalError()); + + for (unsigned int id=0; id!=dim; ++id) + { + error_indicators_vector[present_cell](id) + -= 0.5*face_integrals[cell->face(face_no)](id); + } + + } + + estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm(); + + } + + } + + + // @sect4{PlasticityContactProblem::refine_grid} + + // If you've made it this far into the deal.II tutorial, the following + // function refining the mesh should not pose any challenges to you + // any more. It refines the mesh, either globally or using the Kelly + // error estimator, and if so asked also transfers the solution from + // the previous to the next mesh. In the latter case, we also need + // to compute the active set and other quantities again, for which we + // need the information computed by compute_nonlinear_residual(). + template + void + ElastoPlasticProblem::refine_grid () + { + // --------------------------------------------------------------- + // Make a field variable for history varibales to be able to + // transfer the data to the quadrature points of the new mesh + FE_DGQ history_fe (1); + DoFHandler history_dof_handler (triangulation); + history_dof_handler.distribute_dofs (history_fe); + std::vector< std::vector< Vector > > + history_stress_field (dim, std::vector< Vector >(dim)), + local_history_stress_values_at_qpoints (dim, std::vector< Vector >(dim)), + local_history_stress_fe_values (dim, std::vector< Vector >(dim)); + + + std::vector< std::vector< Vector > > + history_strain_field (dim, std::vector< Vector >(dim)), + local_history_strain_values_at_qpoints (dim, std::vector< Vector >(dim)), + local_history_strain_fe_values (dim, std::vector< Vector >(dim)); + + for (unsigned int i=0; i qpoint_to_dof_matrix (history_fe.dofs_per_cell, + quadrature_formula.size()); + FETools::compute_projection_from_quadrature_points_matrix + (history_fe, + quadrature_formula, quadrature_formula, + qpoint_to_dof_matrix); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + dg_cell = history_dof_handler.begin_active(); + for (; cell!=endc; ++cell, ++dg_cell) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + for (unsigned int i=0; iset_dof_values (local_history_stress_fe_values[i][j], + history_stress_field[i][j]); + + qpoint_to_dof_matrix.vmult (local_history_strain_fe_values[i][j], + local_history_strain_values_at_qpoints[i][j]); + dg_cell->set_dof_values (local_history_strain_fe_values[i][j], + history_strain_field[i][j]); + } + } + + + // --------------------------------------------------------------- + // Refine the mesh + if (refinement_strategy == RefinementStrategy::refine_global) + { + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned()) + cell->set_refine_flag (); + } + else + { + const double refine_fraction_cells = .3, + coarsen_fraction_cells = .03; +// const double refine_fraction_cells = .1, +// coarsen_fraction_cells = .3; + + parallel::distributed::GridRefinement + ::refine_and_coarsen_fixed_number(triangulation, + estimated_error_per_cell, + refine_fraction_cells, coarsen_fraction_cells); + } + + triangulation.prepare_coarsening_and_refinement(); + + parallel::distributed::SolutionTransfer solution_transfer(dof_handler); + solution_transfer.prepare_for_coarsening_and_refinement(solution); + + + parallel::distributed::SolutionTransfer incremental_displacement_transfer(dof_handler); + if (transfer_solution) + incremental_displacement_transfer.prepare_for_coarsening_and_refinement(incremental_displacement); + + SolutionTransfer > history_stress_field_transfer0(history_dof_handler), + history_stress_field_transfer1(history_dof_handler), + history_stress_field_transfer2(history_dof_handler); + history_stress_field_transfer0.prepare_for_coarsening_and_refinement(history_stress_field[0]); + if ( dim > 1) + { + history_stress_field_transfer1.prepare_for_coarsening_and_refinement(history_stress_field[1]); + } + if ( dim == 3) + { + history_stress_field_transfer2.prepare_for_coarsening_and_refinement(history_stress_field[2]); + } + + SolutionTransfer > history_strain_field_transfer0(history_dof_handler), + history_strain_field_transfer1(history_dof_handler), + history_strain_field_transfer2(history_dof_handler); + history_strain_field_transfer0.prepare_for_coarsening_and_refinement(history_strain_field[0]); + if ( dim > 1) + { + history_strain_field_transfer1.prepare_for_coarsening_and_refinement(history_strain_field[1]); + } + if ( dim == 3) + { + history_strain_field_transfer2.prepare_for_coarsening_and_refinement(history_strain_field[2]); + } + + triangulation.execute_coarsening_and_refinement(); + pcout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system(); + setup_quadrature_point_history (); + + + TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs, mpi_communicator); +// distributed_solution = solution; + solution_transfer.interpolate(distributed_solution); + solution = distributed_solution; + + if (transfer_solution) + { + TrilinosWrappers::MPI::Vector distributed_incremental_displacement(locally_owned_dofs, mpi_communicator); +// distributed_incremental_displacement = incremental_displacement; + incremental_displacement_transfer.interpolate(distributed_incremental_displacement); + incremental_displacement = distributed_incremental_displacement; +// compute_nonlinear_residual(incremental_displacement); + } + + // --------------------------------------------------- + history_dof_handler.distribute_dofs (history_fe); + // stress + std::vector< std::vector< Vector > > + distributed_history_stress_field (dim, std::vector< Vector >(dim)); + for (unsigned int i=0; i 1) + { + history_stress_field_transfer1.interpolate(history_stress_field[1], distributed_history_stress_field[1]); + } + if ( dim == 3) + { + history_stress_field_transfer2.interpolate(history_stress_field[2], distributed_history_stress_field[2]); + } + + history_stress_field = distributed_history_stress_field; + + // strain + std::vector< std::vector< Vector > > + distributed_history_strain_field (dim, std::vector< Vector >(dim)); + for (unsigned int i=0; i 1) + { + history_strain_field_transfer1.interpolate(history_strain_field[1], distributed_history_strain_field[1]); + } + if ( dim == 3) + { + history_strain_field_transfer2.interpolate(history_strain_field[2], distributed_history_strain_field[2]); + } + + history_strain_field = distributed_history_strain_field; + + // --------------------------------------------------------------- + // Transfer the history data to the quadrature points of the new mesh + // In a final step, we have to get the data back from the now + // interpolated global field to the quadrature points on the + // new mesh. The following code will do that: + + FullMatrix dof_to_qpoint_matrix (quadrature_formula.size(), + history_fe.dofs_per_cell); + FETools::compute_interpolation_to_quadrature_points_matrix + (history_fe, + quadrature_formula, + dof_to_qpoint_matrix); + cell = dof_handler.begin_active(); + endc = dof_handler.end(); + dg_cell = history_dof_handler.begin_active(); + for (; cell != endc; ++cell, ++dg_cell) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + for (unsigned int i=0; iget_dof_values (history_stress_field[i][j], + local_history_stress_fe_values[i][j]); + dof_to_qpoint_matrix.vmult (local_history_stress_values_at_qpoints[i][j], + local_history_stress_fe_values[i][j]); + + dg_cell->get_dof_values (history_strain_field[i][j], + local_history_strain_fe_values[i][j]); + dof_to_qpoint_matrix.vmult (local_history_strain_values_at_qpoints[i][j], + local_history_strain_fe_values[i][j]); + for (unsigned int q=0; quser_pointer for this that is available in each cell. + // + // To put this into larger perspective, we note that if we had previously + // available stresses in our model (which we assume do not exist for the + // purpose of this program), then we would need to interpolate the field of + // preexisting stresses to the quadrature points. Likewise, if we were to + // simulate elasto-plastic materials with hardening/softening, then we would + // have to store additional history variables like the present yield stress + // of the accumulated plastic strains in each quadrature + // points. Pre-existing hardening or weakening would then be implemented by + // interpolating these variables in the present function as well. + template + void ElastoPlasticProblem::setup_quadrature_point_history () + { + // What we need to do here is to first count how many quadrature points + // are within the responsibility of this processor. This, of course, + // equals the number of cells that belong to this processor times the + // number of quadrature points our quadrature formula has on each cell. + // + // For good measure, we also set all user pointers of all cells, whether + // ours of not, to the null pointer. This way, if we ever access the user + // pointer of a cell which we should not have accessed, a segmentation + // fault will let us know that this should not have happened: + unsigned int our_cells = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned()) + ++our_cells; + + triangulation.clear_user_data(); + + // Next, allocate as many quadrature objects as we need. Since the + // resize function does not actually shrink the amount of + // allocated memory if the requested new size is smaller than the old + // size, we resort to a trick to first free all memory, and then + // reallocate it: we declare an empty vector as a temporary variable and + // then swap the contents of the old vector and this temporary + // variable. This makes sure that the + // quadrature_point_history is now really empty, and we can + // let the temporary variable that now holds the previous contents of the + // vector go out of scope and be destroyed. In the next step. we can then + // re-allocate as many elements as we need, with the vector + // default-initializing the PointHistory objects, which + // includes setting the stress variables to zero. + { + std::vector > tmp; + tmp.swap (quadrature_point_history); + } + quadrature_point_history.resize (our_cells * + quadrature_formula.size()); + + // Finally loop over all cells again and set the user pointers from the + // cells that belong to the present processor to point to the first + // quadrature point objects corresponding to this cell in the vector of + // such objects: + unsigned int history_index = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned()) + { + cell->set_user_pointer (&quadrature_point_history[history_index]); + history_index += quadrature_formula.size(); + } + + // At the end, for good measure make sure that our count of elements was + // correct and that we have both used up all objects we allocated + // previously, and not point to any objects beyond the end of the + // vector. Such defensive programming strategies are always good checks to + // avoid accidental errors and to guard against future changes to this + // function that forget to update all uses of a variable at the same + // time. Recall that constructs using the Assert macro are + // optimized away in optimized mode, so do not affect the run time of + // optimized runs: + Assert (history_index == quadrature_point_history.size(), + ExcInternalError()); + } + + // @sect4{ElastoPlasticProblem::update_quadrature_point_history} + + // At the end of each time step, we should have computed an incremental + // displacement update so that the material in its new configuration + // accommodates for the difference between the external body and boundary + // forces applied during this time step minus the forces exerted through + // preexisting internal stresses. In order to have the preexisting + // stresses available at the next time step, we therefore have to update the + // preexisting stresses with the stresses due to the incremental + // displacement computed during the present time step. Ideally, the + // resulting sum of internal stresses would exactly counter all external + // forces. Indeed, a simple experiment can make sure that this is so: if we + // choose boundary conditions and body forces to be time independent, then + // the forcing terms (the sum of external forces and internal stresses) + // should be exactly zero. If you make this experiment, you will realize + // from the output of the norm of the right hand side in each time step that + // this is almost the case: it is not exactly zero, since in the first time + // step the incremental displacement and stress updates were computed + // relative to the undeformed mesh, which was then deformed. In the second + // time step, we again compute displacement and stress updates, but this + // time in the deformed mesh -- there, the resulting updates are very small + // but not quite zero. This can be iterated, and in each such iteration the + // residual, i.e. the norm of the right hand side vector, is reduced; if one + // makes this little experiment, one realizes that the norm of this residual + // decays exponentially with the number of iterations, and after an initial + // very rapid decline is reduced by roughly a factor of about 3.5 in each + // iteration (for one testcase I looked at, other testcases, and other + // numbers of unknowns change the factor, but not the exponential decay). + + // In a sense, this can then be considered as a quasi-timestepping scheme to + // resolve the nonlinear problem of solving large-deformation elasticity on + // a mesh that is moved along in a Lagrangian manner. + // + // Another complication is that the existing (old) stresses are defined on + // the old mesh, which we will move around after updating the stresses. If + // this mesh update involves rotations of the cell, then we need to also + // rotate the updated stress, since it was computed relative to the + // coordinate system of the old cell. + // + // Thus, what we need is the following: on each cell which the present + // processor owns, we need to extract the old stress from the data stored + // with each quadrature point, compute the stress update, add the two + // together, and then rotate the result together with the incremental + // rotation computed from the incremental displacement at the present + // quadrature point. We will detail these steps below: + template + void ElastoPlasticProblem:: + update_quadrature_point_history () + { + // First, set up an FEValues object by which we will evaluate + // the displacements and the gradients thereof at the + // quadrature points, together with a vector that will hold this + // information: + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points); + + const unsigned int n_q_points = quadrature_formula.size(); + + std::vector > incremental_strain_tensor(n_q_points); + SymmetricTensor<4, dim> stress_strain_tensor; + + + // Then loop over all cells and do the job in the cells that belong to our + // subdomain: + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell != endc; ++cell) + if (cell->is_locally_owned()) + { + // Next, get a pointer to the quadrature point history data local to + // the present cell, and, as a defensive measure, make sure that + // this pointer is within the bounds of the global array: + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + + // Then initialize the FEValues object on the present + // cell, and extract the strains of the displacement at the + // quadrature points + fe_values.reinit (cell); + fe_values[displacement].get_function_symmetric_gradients(incremental_displacement, + incremental_strain_tensor); + + // Then loop over the quadrature points of this cell: + for (unsigned int q=0; qrun() + // have to do with generating output. The following one is an attempt + // at showing the deformed body in its deformed configuration. To this + // end, this function takes a displacement vector field and moves every + // vertex of the (local part) of the mesh by the previously computed + // displacement. We will call this function with the current + // displacement field before we generate graphical output, and we will + // call it again after generating graphical output with the negative + // displacement field to undo the changes to the mesh so made. + // + // The function itself is pretty straightforward. All we have to do + // is keep track which vertices we have already touched, as we + // encounter the same vertices multiple times as we loop over cells. + template + void + ElastoPlasticProblem:: + move_mesh (const TrilinosWrappers::MPI::Vector &displacement) const + { + std::vector vertex_touched(triangulation.n_vertices(), false); + + for (typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned()) + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + if (vertex_touched[cell->vertex_index(v)] == false) + { + vertex_touched[cell->vertex_index(v)] = true; + + Point vertex_displacement; + for (unsigned int d = 0; d < dim; ++d) + vertex_displacement[d] = displacement(cell->vertex_dof_index(v, d)); + + cell->vertex(v) += vertex_displacement; + } + } + + + + // @sect4{PlasticityContactProblem::output_results} + + // Next is the function we use to actually generate graphical output. The + // function is a bit tedious, but not actually particularly complicated. + // It moves the mesh at the top (and moves it back at the end), then + // computes the contact forces along the contact surface. We can do + // so (as shown in the accompanying paper) by taking the untreated + // residual vector and identifying which degrees of freedom + // correspond to those with contact by asking whether they have an + // inhomogeneous constraints associated with them. As always, we need + // to be mindful that we can only write into completely distributed + // vectors (i.e., vectors without ghost elements) but that when we + // want to generate output, we need vectors that do indeed have + // ghost entries for all locally relevant degrees of freedom. + template + void + ElastoPlasticProblem::output_results (const std::string &filename_base) + { + TimerOutput::Scope t(computing_timer, "Graphical output"); + + pcout << " Writing graphical output... " << std::flush; + + TrilinosWrappers::MPI::Vector magnified_solution(solution); + + const double magnified_factor = 3; + magnified_solution *= magnified_factor; + + move_mesh(magnified_solution); + + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + + // + const std::vector + data_component_interpretation(dim, DataComponentInterpretation::component_is_part_of_vector); + data_out.add_data_vector(solution, + std::vector (dim, "displacement"), + DataOut::type_dof_data, data_component_interpretation); + + // + std::vector solution_names; + + switch (dim) + { + case 1: + solution_names.push_back ("displacement"); + break; + case 2: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + break; + case 3: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + solution_names.push_back ("z_displacement"); + break; + default: + AssertThrow (false, ExcNotImplemented()); + } + + data_out.add_data_vector (solution, solution_names); + + + // + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = triangulation.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + // + data_out.add_data_vector(fraction_of_plastic_q_points_per_cell, + "fraction_of_plastic_q_points"); + + // + data_out.build_patches(); + + // In the remainder of the function, we generate one VTU file on + // every processor, indexed by the subdomain id of this processor. + // On the first processor, we then also create a .pvtu + // file that indexes all of the VTU files so that the entire + // set of output files can be read at once. These .pvtu + // are used by Paraview to describe an entire parallel computation's + // output files. We then do the same again for the competitor of + // Paraview, the Visit visualization program, by creating a matching + // .visit file. + const std::string filename = + (output_dir + filename_base + "-" + + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4)); + + std::ofstream output_vtu((filename + ".vtu").c_str()); + data_out.write_vtu(output_vtu); + pcout << output_dir + filename_base << ".pvtu" << std::endl; + + + if (this_mpi_process == 0) + { + std::vector filenames; + for (unsigned int i = 0; i < n_mpi_processes; ++i) + filenames.push_back(filename_base + "-" + + Utilities::int_to_string(i, 4) + + ".vtu"); + + std::ofstream pvtu_master_output((output_dir + filename_base + ".pvtu").c_str()); + data_out.write_pvtu_record(pvtu_master_output, filenames); + + std::ofstream visit_master_output((output_dir + filename_base + ".visit").c_str()); + data_out.write_visit_record(visit_master_output, filenames); + + // produce eps files for mesh illustration + std::ofstream output_eps((filename + ".eps").c_str()); + GridOut grid_out; + grid_out.write_eps(triangulation, output_eps); + } + + // Extrapolate the stresses from Gauss point to the nodes + SymmetricTensor<2, dim> stress_at_qpoint; + + FE_DGQ history_fe (1); + DoFHandler history_dof_handler (triangulation); + history_dof_handler.distribute_dofs (history_fe); + std::vector< std::vector< Vector > > + history_stress_field (dim, std::vector< Vector >(dim)), + local_history_stress_values_at_qpoints (dim, std::vector< Vector >(dim)), + local_history_stress_fe_values (dim, std::vector< Vector >(dim)); + for (unsigned int i=0; i VM_stress_field (history_dof_handler.n_dofs()), + local_VM_stress_values_at_qpoints (quadrature_formula.size()), + local_VM_stress_fe_values (history_fe.dofs_per_cell); + + FullMatrix qpoint_to_dof_matrix (history_fe.dofs_per_cell, + quadrature_formula.size()); + FETools::compute_projection_from_quadrature_points_matrix + (history_fe, + quadrature_formula, quadrature_formula, + qpoint_to_dof_matrix); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + dg_cell = history_dof_handler.begin_active(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell!=endc; ++cell, ++dg_cell) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + + // Then loop over the quadrature points of this cell: + for (unsigned int q=0; qset_dof_values (local_history_stress_fe_values[i][j], + history_stress_field[i][j]); + } + + qpoint_to_dof_matrix.vmult (local_VM_stress_fe_values, + local_VM_stress_values_at_qpoints); + dg_cell->set_dof_values (local_VM_stress_fe_values, + VM_stress_field); + + + } + + // Save stresses on nodes by nodal averaging + // construct a DoFHandler object based on FE_Q with 1 degree of freedom + // in order to compute stresses on nodes (by applying nodal averaging) + // Therefore, each vertex has one degree of freedom + FE_Q fe_1 (1); + DoFHandler dof_handler_1 (triangulation); + dof_handler_1.distribute_dofs (fe_1); + + AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(), + ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices())); + + std::vector< std::vector< Vector > > + history_stress_on_vertices (dim, std::vector< Vector >(dim)); + for (unsigned int i=0; i VM_stress_on_vertices (dof_handler_1.n_dofs()), + counter_on_vertices (dof_handler_1.n_dofs()); + VM_stress_on_vertices = 0; + counter_on_vertices = 0; + + cell = dof_handler.begin_active(); + dg_cell = history_dof_handler.begin_active(); + typename DoFHandler::active_cell_iterator + cell_1 = dof_handler_1.begin_active(); + for (; cell!=endc; ++cell, ++dg_cell, ++cell_1) + if (cell->is_locally_owned()) + { + dg_cell->get_dof_values (VM_stress_field, + local_VM_stress_fe_values); + + for (unsigned int i=0; iget_dof_values (history_stress_field[i][j], + local_history_stress_fe_values[i][j]); + } + + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + { + types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0); + + // begin check + // Point point1, point2; + // point1 = cell_1->vertex(v); + // point2 = dg_cell->vertex(v); + // AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError()); + // end check + + counter_on_vertices (dof_1_vertex) += 1; + + VM_stress_on_vertices (dof_1_vertex) += local_VM_stress_fe_values (v); + + for (unsigned int i=0; i data_out; + data_out.attach_dof_handler (history_dof_handler); + + + data_out.add_data_vector (history_stress_field[0][0], "stress_xx"); + data_out.add_data_vector (history_stress_field[1][1], "stress_yy"); + data_out.add_data_vector (history_stress_field[0][1], "stress_xy"); + data_out.add_data_vector (VM_stress_field, "Von_Mises_stress"); + + if (dim == 3) + { + data_out.add_data_vector (history_stress_field[0][2], "stress_xz"); + data_out.add_data_vector (history_stress_field[1][2], "stress_yz"); + data_out.add_data_vector (history_stress_field[2][2], "stress_zz"); + } + + data_out.build_patches (); + + const std::string filename_base_stress = ("stress-" + filename_base); + + const std::string filename = + (output_dir + filename_base_stress + "-" + + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4)); + + std::ofstream output_vtu((filename + ".vtu").c_str()); + data_out.write_vtu(output_vtu); + pcout << output_dir + filename_base_stress << ".pvtu" << std::endl; + + if (this_mpi_process == 0) + { + std::vector filenames; + for (unsigned int i = 0; i < n_mpi_processes; ++i) + filenames.push_back(filename_base_stress + "-" + + Utilities::int_to_string(i, 4) + + ".vtu"); + + std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str()); + data_out.write_pvtu_record(pvtu_master_output, filenames); + + std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str()); + data_out.write_visit_record(visit_master_output, filenames); + } + + + } + + { + DataOut data_out; + data_out.attach_dof_handler (dof_handler_1); + + + data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_xx_averaged"); + data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_yy_averaged"); + data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_xy_averaged"); + data_out.add_data_vector (VM_stress_on_vertices, "Von_Mises_stress_averaged"); + + if (dim == 3) + { + data_out.add_data_vector (history_stress_on_vertices[0][2], "stress_xz_averaged"); + data_out.add_data_vector (history_stress_on_vertices[1][2], "stress_yz_averaged"); + data_out.add_data_vector (history_stress_on_vertices[2][2], "stress_zz_averaged"); + } + + data_out.build_patches (); + + const std::string filename_base_stress = ("averaged-stress-" + filename_base); + + const std::string filename = + (output_dir + filename_base_stress + "-" + + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4)); + + std::ofstream output_vtu((filename + ".vtu").c_str()); + data_out.write_vtu(output_vtu); + pcout << output_dir + filename_base_stress << ".pvtu" << std::endl; + + if (this_mpi_process == 0) + { + std::vector filenames; + for (unsigned int i = 0; i < n_mpi_processes; ++i) + filenames.push_back(filename_base_stress + "-" + + Utilities::int_to_string(i, 4) + + ".vtu"); + + std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str()); + data_out.write_pvtu_record(pvtu_master_output, filenames); + + std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str()); + data_out.write_visit_record(visit_master_output, filenames); + } + + + } + // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + } + + magnified_solution *= -1; + move_mesh(magnified_solution); + + // Timoshenko beam + if (base_mesh == "Timoshenko beam") + { + const double length = .48, + depth = .12; + + Point intersted_point(length, -depth/2); + Point vertex_displacement; + bool vertex_found = false; + + for (typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned() && !vertex_found) + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + if ( std::fabs(cell->vertex(v)[0] - intersted_point[0])<1e-6 && + std::fabs(cell->vertex(v)[1] - intersted_point[1])<1e-6) + { + vertex_found = true; + + for (unsigned int d = 0; d < dim; ++d) + vertex_displacement[d] = solution(cell->vertex_dof_index(v, d)); + + break; + } + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + AssertThrow(vertex_found, ExcInternalError()); + std::cout << "Displacement at the point (" << intersted_point[0] + << ", " << intersted_point[1] << ") is " + << "(" << vertex_displacement[0] + << ", " << vertex_displacement[1] << ").\n"; + + Vector vertex_exact_displacement(dim); + EquationData::IncrementalBoundaryValues incremental_boundary_values(present_time, end_time); + incremental_boundary_values.vector_value (intersted_point, vertex_exact_displacement); + + std::cout << "Exact displacement at the point (" << intersted_point[0] + << ", " << intersted_point[1] << ") is " + << "(" << vertex_exact_displacement[0] + << ", " << vertex_exact_displacement[1] << ").\n\n"; + + }else if (base_mesh == "Thick_tube_internal_pressure") + { + const double pressure (0.6*2.4e8), + inner_radius (.1); +// const double pressure (1.94e8), +// inner_radius (.1); + + + // Plane stress +// const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))); + // 3d and plane strain + const double mu (e_modulus / (2 * (1 + nu))); + + const Point point_A(inner_radius, 0.); + Vector disp_A(dim); + + // make a non-parallel copy of solution + Vector copy_solution(solution); + + typename Evaluation::PointValuesEvaluation:: + PointValuesEvaluation point_values_evaluation(point_A); + + point_values_evaluation.compute (dof_handler, copy_solution, disp_A); + + table_results.add_value("time step", timestep_no); + table_results.add_value("Cells", triangulation.n_global_active_cells()); + table_results.add_value("DoFs", dof_handler.n_dofs()); + table_results.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0); + table_results.add_value("4*mu*u_A/(sigma_0*a)", 4*mu*disp_A(0)/(sigma_0*inner_radius)); + + // Compute stresses in the POLAR coordinates, 1- save it on Gauss points, + // 2- extrapolate them to nodes and taking their avarages (nodal avaraging) + AssertThrow (dim == 2, ExcNotImplemented()); + + // we define a rotation matrix to be able to transform the stress + // from the Cartesian coordinate to the polar coordinate + Tensor<2, dim> rotation_matrix; // [cos sin; -sin cos] , sigma_r = rot * sigma * rot^T + + FEValues fe_values (fe, quadrature_formula, update_quadrature_points | + update_values | update_gradients); + + const unsigned int n_q_points = quadrature_formula.size(); + + std::vector > strain_tensor(n_q_points); + SymmetricTensor<4, dim> stress_strain_tensor; + Tensor<2, dim> stress_at_qpoint; + + FE_DGQ history_fe (1); + DoFHandler history_dof_handler (triangulation); + history_dof_handler.distribute_dofs (history_fe); + std::vector< std::vector< Vector > > + history_stress_field (dim, std::vector< Vector >(dim)), + local_history_stress_values_at_qpoints (dim, std::vector< Vector >(dim)), + local_history_stress_fe_values (dim, std::vector< Vector >(dim)); + for (unsigned int i=0; i qpoint_to_dof_matrix (history_fe.dofs_per_cell, + quadrature_formula.size()); + FETools::compute_projection_from_quadrature_points_matrix + (history_fe, + quadrature_formula, quadrature_formula, + qpoint_to_dof_matrix); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + dg_cell = history_dof_handler.begin_active(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell!=endc; ++cell, ++dg_cell) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + + // Then loop over the quadrature points of this cell: + for (unsigned int q=0; q point = local_quadrature_points_history[q].point; + const double radius = point.norm (); + const double theta = std::atan2(point(1),point(0)); + + // rotation matrix + rotation_matrix[0][0] = std::cos(theta); + rotation_matrix[0][1] = std::sin(theta); + rotation_matrix[1][0] = -std::sin(theta); + rotation_matrix[1][1] = std::cos(theta); + + // stress in polar coordinate + stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix); + + for (unsigned int i=0; iset_dof_values (local_history_stress_fe_values[i][j], + history_stress_field[i][j]); + } + + } + + { + DataOut data_out; + data_out.attach_dof_handler (history_dof_handler); + + + data_out.add_data_vector (history_stress_field[0][0], "stress_rr"); + data_out.add_data_vector (history_stress_field[1][1], "stress_tt"); + data_out.add_data_vector (history_stress_field[0][1], "stress_rt"); + + data_out.build_patches (); + + const std::string filename_base_stress = ("stress-polar-" + filename_base); + + const std::string filename = + (output_dir + filename_base_stress + "-" + + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4)); + + std::ofstream output_vtu((filename + ".vtu").c_str()); + data_out.write_vtu(output_vtu); + pcout << output_dir + filename_base_stress << ".pvtu" << std::endl; + + if (this_mpi_process == 0) + { + std::vector filenames; + for (unsigned int i = 0; i < n_mpi_processes; ++i) + filenames.push_back(filename_base_stress + "-" + + Utilities::int_to_string(i, 4) + + ".vtu"); + + std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str()); + data_out.write_pvtu_record(pvtu_master_output, filenames); + + std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str()); + data_out.write_visit_record(visit_master_output, filenames); + } + + + } + + // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + // construct a DoFHandler object based on FE_Q with 1 degree of freedom + // in order to compute stresses on nodes (by applying nodal averaging) + // Therefore, each vertex has one degree of freedom + FE_Q fe_1 (1); + DoFHandler dof_handler_1 (triangulation); + dof_handler_1.distribute_dofs (fe_1); + + AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(), + ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices())); + + std::vector< std::vector< Vector > > + history_stress_on_vertices (dim, std::vector< Vector >(dim)); + for (unsigned int i=0; i counter_on_vertices (dof_handler_1.n_dofs()); + counter_on_vertices = 0; + + cell = dof_handler.begin_active(); + dg_cell = history_dof_handler.begin_active(); + typename DoFHandler::active_cell_iterator + cell_1 = dof_handler_1.begin_active(); + for (; cell!=endc; ++cell, ++dg_cell, ++cell_1) + if (cell->is_locally_owned()) + { + + for (unsigned int i=0; iget_dof_values (history_stress_field[i][j], + local_history_stress_fe_values[i][j]); + } + + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + { + types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0); + + // begin check +// Point point1, point2; +// point1 = cell_1->vertex(v); +// point2 = dg_cell->vertex(v); +// AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError()); + // end check + + counter_on_vertices (dof_1_vertex) += 1; + + for (unsigned int i=0; i data_out; + data_out.attach_dof_handler (dof_handler_1); + + + data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_rr_averaged"); + data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_tt_averaged"); + data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_rt_averaged"); + + data_out.build_patches (); + + const std::string filename_base_stress = ("averaged-stress-polar-" + filename_base); + + const std::string filename = + (output_dir + filename_base_stress + "-" + + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4)); + + std::ofstream output_vtu((filename + ".vtu").c_str()); + data_out.write_vtu(output_vtu); + pcout << output_dir + filename_base_stress << ".pvtu" << std::endl; + + if (this_mpi_process == 0) + { + std::vector filenames; + for (unsigned int i = 0; i < n_mpi_processes; ++i) + filenames.push_back(filename_base_stress + "-" + + Utilities::int_to_string(i, 4) + + ".vtu"); + + std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str()); + data_out.write_pvtu_record(pvtu_master_output, filenames); + + std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str()); + data_out.write_visit_record(visit_master_output, filenames); + } + + + } + // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + if ( std::abs( (present_time/end_time)*(pressure/sigma_0) - 0.6 ) < + .501*(present_timestep/end_time)*(pressure/sigma_0) ) + { + + // table_results_2: presenting the stress_rr and stress_tt on the nodes of bottom edge + const unsigned int face_id = 3; + + std::vector vertices_found (dof_handler_1.n_dofs(), false); + + bool evaluation_face_found = false; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + cell_1 = dof_handler_1.begin_active(); + for (; cell!=endc; ++cell, ++cell_1) + if (cell->is_locally_owned()) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && + cell->face(face)->boundary_indicator() == face_id) + { + if (!evaluation_face_found) + { + evaluation_face_found = true; + } + + + for (unsigned int v = 0; v < GeometryInfo::vertices_per_face; ++v) + { + types::global_dof_index dof_1_vertex = + cell_1->face(face)->vertex_dof_index(v, 0); + if (!vertices_found[dof_1_vertex]) + { + + const Point vertex_coordinate = cell_1->face(face)->vertex(v); + + table_results_2.add_value("x coordinate", vertex_coordinate[0]); + table_results_2.add_value("stress_rr", history_stress_on_vertices[0][0](dof_1_vertex)); + table_results_2.add_value("stress_tt", history_stress_on_vertices[1][1](dof_1_vertex)); + table_results_2.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0); + + vertices_found[dof_1_vertex] = true; + } + } + + } + } + + } + + AssertThrow(evaluation_face_found, ExcInternalError()); + + // table_results_3: presenting the mean stress_rr of the nodes on the inner radius + const unsigned int face_id_2 = 0; + + Tensor<2, dim> stress_node, + mean_stress_polar; + mean_stress_polar = 0; + + std::vector vertices_found_2 (dof_handler_1.n_dofs(), false); + unsigned int no_vertices_found = 0; + + evaluation_face_found = false; + + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + cell_1 = dof_handler_1.begin_active(); + for (; cell!=endc; ++cell, ++cell_1) + if (cell->is_locally_owned()) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && + cell->face(face)->boundary_indicator() == face_id_2) + { + if (!evaluation_face_found) + { + evaluation_face_found = true; + } + + + for (unsigned int v = 0; v < GeometryInfo::vertices_per_face; ++v) + { + types::global_dof_index dof_1_vertex = + cell_1->face(face)->vertex_dof_index(v, 0); + if (!vertices_found_2[dof_1_vertex]) + { + for (unsigned int ir=0; ir point_C(0., height); + Vector disp_C(dim); + + // make a non-parallel copy of solution + Vector copy_solution(solution); + + typename Evaluation::PointValuesEvaluation:: + PointValuesEvaluation point_values_evaluation(point_C); + + point_values_evaluation.compute (dof_handler, copy_solution, disp_C); + + table_results.add_value("time step", timestep_no); + table_results.add_value("Cells", triangulation.n_global_active_cells()); + table_results.add_value("DoFs", dof_handler.n_dofs()); + table_results.add_value("4*mu*u_C/(sigma_0*r)", 4*mu*disp_C(1)/(sigma_0*inner_radius)); + } + */ + + // compute average sigma_yy on the bottom edge + double stress_yy_av; + { + stress_yy_av = 0; + const unsigned int face_id = 1; + + std::vector vertices_found (dof_handler_1.n_dofs(), false); + unsigned int no_vertices_in_face = 0; + + bool evaluation_face_found = false; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + cell_1 = dof_handler_1.begin_active(); + for (; cell!=endc; ++cell, ++cell_1) + if (cell->is_locally_owned()) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && + cell->face(face)->boundary_indicator() == face_id) + { + if (!evaluation_face_found) + { + evaluation_face_found = true; + } + + + for (unsigned int v = 0; v < GeometryInfo::vertices_per_face; ++v) + { + types::global_dof_index dof_1_vertex = + cell_1->face(face)->vertex_dof_index(v, 0); + if (!vertices_found[dof_1_vertex]) + { + stress_yy_av += history_stress_on_vertices[1][1](dof_1_vertex); + ++no_vertices_in_face; + + vertices_found[dof_1_vertex] = true; + } + } + + } + } + + } + + AssertThrow(evaluation_face_found, ExcInternalError()); + + stress_yy_av /= no_vertices_in_face; + + } + + // table_results_2: Demonstrate the stress_yy on the nodes of bottom edge + +// if ( std::abs( (stress_yy_av/sigma_0) - .91 ) < .2 ) + if ( (timestep_no) % 19 == 0 ) +// if ( true ) + { + const unsigned int face_id = 1; + + std::vector vertices_found (dof_handler_1.n_dofs(), false); + + bool evaluation_face_found = false; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + cell_1 = dof_handler_1.begin_active(); + for (; cell!=endc; ++cell, ++cell_1) + if (cell->is_locally_owned()) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && + cell->face(face)->boundary_indicator() == face_id) + { + if (!evaluation_face_found) + { + evaluation_face_found = true; + } + + + for (unsigned int v = 0; v < GeometryInfo::vertices_per_face; ++v) + { + types::global_dof_index dof_1_vertex = + cell_1->face(face)->vertex_dof_index(v, 0); + + const Point vertex_coordinate = cell_1->face(face)->vertex(v); + + if (!vertices_found[dof_1_vertex] && std::abs(vertex_coordinate[2])<1.e-8) + { + table_results_2.add_value("x", vertex_coordinate[0]); + table_results_2.add_value("x/r", vertex_coordinate[0]/inner_radius); + table_results_2.add_value("stress_xx/sigma_0", history_stress_on_vertices[0][0](dof_1_vertex)/sigma_0); + table_results_2.add_value("stress_yy/sigma_0", history_stress_on_vertices[1][1](dof_1_vertex)/sigma_0); + table_results_2.add_value("stress_yy_av/sigma_0", stress_yy_av/sigma_0); + table_results_2.add_value("Imposed u_y", (imposed_displacement*present_time/end_time)); + + vertices_found[dof_1_vertex] = true; + } + } + + } + } + + } + + AssertThrow(evaluation_face_found, ExcInternalError()); + + } + + // table_results_3: Demonstrate the Stress_mean (average tensile stress) + // on the bottom edge versus epsilon_yy on the bottom left corner + { + double strain_yy_A; + + // compute strain_yy_A + // Since the point A is the node on the bottom left corner, + // we need to work just with one element + { + const Point point_A(inner_radius, 0, 0); + + Vector local_strain_yy_values_at_qpoints (quadrature_formula.size()), + local_strain_yy_fe_values (history_fe.dofs_per_cell); + + SymmetricTensor<2, dim> strain_at_qpoint; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + dg_cell = history_dof_handler.begin_active(); + + bool cell_found = false; + + for (; cell!=endc; ++cell, ++dg_cell) + if (cell->is_locally_owned() && !cell_found) + { + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + if ( std::fabs(cell->vertex(v)[0] - point_A[0])<1e-6 && + std::fabs(cell->vertex(v)[1] - point_A[1])<1e-6 && + std::fabs(cell->vertex(v)[2] - point_A[2])<1e-6) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + + // Then loop over the quadrature points of this cell: + for (unsigned int q=0; q point_A(0, height/2, length); + Vector disp_A(dim); + + // make a non-parallel copy of solution + Vector copy_solution(solution); + + typename Evaluation::PointValuesEvaluation:: + PointValuesEvaluation point_values_evaluation(point_A); + + point_values_evaluation.compute (dof_handler, copy_solution, disp_A); + + table_results.add_value("time step", timestep_no); + table_results.add_value("Cells", triangulation.n_global_active_cells()); + table_results.add_value("DoFs", dof_handler.n_dofs()); + table_results.add_value("pressure", pressure*present_time/end_time); + table_results.add_value("u_A", disp_A(1)); + } + + { + // demonstrate the location and maximum von-Mises stress in the + // specified domain close to the clamped face, z = 0 + // top domain: height/2 - thickness_flange <= y <= height/2 + // 0 <= z <= 2 * thickness_flange + // bottom domain: -height/2 <= y <= -height/2 + thickness_flange + // 0 <= z <= 2 * thickness_flange + + double VM_stress_max (0); + Point point_max; + + SymmetricTensor<2, dim> stress_at_qpoint; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + const FEValuesExtractors::Vector displacement(0); + + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(cell->user_pointer()); + Assert (local_quadrature_points_history >= + &quadrature_point_history.front(), + ExcInternalError()); + Assert (local_quadrature_points_history < + &quadrature_point_history.back(), + ExcInternalError()); + + // Then loop over the quadrature points of this cell: + for (unsigned int q=0; q VM_stress_max) + { + VM_stress_max = VM_stress; + point_max = local_quadrature_points_history[q].point; + } + + } + } + + table_results.add_value("maximum von_Mises stress", VM_stress_max); + table_results.add_value("x", point_max[0]); + table_results.add_value("y", point_max[1]); + table_results.add_value("z", point_max[2]); + + } + + } + + + } + + + // @sect4{PlasticityContactProblem::run} + + // As in all other tutorial programs, the run() function contains + // the overall logic. There is not very much to it here: in essence, it + // performs the loops over all mesh refinement cycles, and within each, hands + // things over to the Newton solver in solve_newton() on the + // current mesh and calls the function that creates graphical output for + // the so-computed solution. It then outputs some statistics concerning both + // run times and memory consumption that has been collected over the course of + // computations on this mesh. + template + void + ElastoPlasticProblem::run () + { + computing_timer.reset(); + + present_time = 0; + present_timestep = 1; + end_time = 10; + timestep_no = 0; + + make_grid(); + + // ---------------------------------------------------------------- + // base_mesh == "Thick_tube_internal_pressure" + /* + const Point center(0, 0); + const double inner_radius = .1, + outer_radius = .2; + + const HyperBallBoundary inner_boundary_description(center, inner_radius); + triangulation.set_boundary (0, inner_boundary_description); + + const HyperBallBoundary outer_boundary_description(center, outer_radius); + triangulation.set_boundary (1, outer_boundary_description); + */ + // ---------------------------------------------------------------- + // base_mesh == "Perforated_strip_tension" + /* + const double inner_radius = 0.05; + + const CylinderBoundary inner_boundary_description(inner_radius, 2); + triangulation.set_boundary (10, inner_boundary_description); + */ + // ---------------------------------------------------------------- + + setup_quadrature_point_history (); + + while (present_time < end_time) + { + present_time += present_timestep; + ++timestep_no; + + if (present_time > end_time) + { + present_timestep -= (present_time - end_time); + present_time = end_time; + } + pcout << std::endl; + pcout << "Time step " << timestep_no << " at time " << present_time + << std::endl; + + relative_error = max_relative_error * 10; + current_refinement_cycle = 0; + + setup_system(); + + + // ------------------------ Refinement based on the relative error ------------------------------- + + while (relative_error >= max_relative_error) + { + solve_newton(); + compute_error(); + + if ( (timestep_no > 1) && (current_refinement_cycle>0) && (relative_error >= max_relative_error) ) + { + pcout << "The relative error, " << relative_error + << " , is still more than maximum relative error, " + << max_relative_error << ", but we move to the next increment.\n"; + relative_error = .1 * max_relative_error; + } + + if (relative_error >= max_relative_error) + { + TimerOutput::Scope t(computing_timer, "Setup: refine mesh"); + ++current_refinement_cycle; + refine_grid(); + } + + } + + // ------------------------ Refinement based on the number of refinement -------------------------- + /* + bool continue_loop = true; + while (continue_loop) + { + solve_newton(); + compute_error(); + + if ( (timestep_no == 1) && (current_refinement_cycle < 1) ) + { + TimerOutput::Scope t(computing_timer, "Setup: refine mesh"); + ++current_refinement_cycle; + refine_grid(); + }else + { + continue_loop = false; + } + + } + */ + + // ------------------------------------------------------------------------------------------------- + + solution += incremental_displacement; + + update_quadrature_point_history (); + + output_results((std::string("solution-") + + Utilities::int_to_string(timestep_no, 4)).c_str()); + + computing_timer.print_summary(); + computing_timer.reset(); + + Utilities::System::MemoryStats stats; + Utilities::System::get_memory_stats(stats); + pcout << "Peak virtual memory used, resident in kB: " << stats.VmSize << " " + << stats.VmRSS << std::endl; + + + if (std::abs(present_time-end_time) < 1.e-7) + { + const std::string filename = (output_dir + "Results"); + + std::ofstream output_txt((filename + ".txt").c_str()); + + pcout << std::endl; + table_results.write_text(output_txt); + pcout << std::endl; + table_results_2.write_text(output_txt); + pcout << std::endl; + table_results_3.write_text(output_txt); + pcout << std::endl; + } + + } + + if (base_mesh == "Thick_tube_internal_pressure") + { + triangulation.set_boundary (0); + triangulation.set_boundary (1); + }else if (base_mesh == "Perforated_strip_tension") + { + triangulation.set_boundary (10); + } + + } +} + +// @sect3{The main function} + +// There really isn't much to the main() function. It looks +// like they always do: +int main (int argc, char *argv[]) +{ + using namespace dealii; + using namespace ElastoPlastic; + + try + { + deallog.depth_console(0); + ParameterHandler prm; + const int dim = 3; + ElastoPlasticProblem::declare_parameters(prm); + if (argc != 2) + { + std::cerr << "*** Call this program as <./elastoplastic input.prm>" << std::endl; + return 1; + } + + prm.read_input(argv[1]); + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + { + ElastoPlasticProblem problem(prm); + problem.run(); + } + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/goal_oriented_elastoplasticity/readme.md b/goal_oriented_elastoplasticity/readme.md new file mode 100644 index 0000000..3b38c22 --- /dev/null +++ b/goal_oriented_elastoplasticity/readme.md @@ -0,0 +1,11 @@ +3d goal-oriented mesh adaptivity in elastoplasticity problems + +The code deals with solving an elastoplasticity problem with linear isotropic hardening. At each load/displacement step, the error based on a prescribed quantity of interest (Goal-oriented error estimation) is computed by using the dual-weighted residual method. + +Based on a prescribed error bound and estimated elementwise errors, the mesh is refined/coarsened. Afterwards, the solution is projected to the new mesh and the analysis process is proceeded. + +The applied methodology and the solved numerical examples can be found in the following paper: + +Ghorashi SSh, Rabczuk T. Goal-Oriented Error Estimation and Mesh Adaptivity in +3d Elastoplasticity Problems. International Journal of Fracture. Accepted. +