From 64cf01451d66906bb78f010d4048269f6fb30083 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 26 May 2016 18:21:13 -0500 Subject: [PATCH] Use the usual indentation style. --- .../elastoplastic.cc | 9509 +++++++++-------- 1 file changed, 4786 insertions(+), 4723 deletions(-) diff --git a/goal_oriented_elastoplasticity/elastoplastic.cc b/goal_oriented_elastoplasticity/elastoplastic.cc index 4b0d9de..449a459 100644 --- a/goal_oriented_elastoplasticity/elastoplastic.cc +++ b/goal_oriented_elastoplasticity/elastoplastic.cc @@ -16,7 +16,7 @@ * * Authors: Seyed Shahram Ghorashi, Bauhaus-Universit\"at Weimar, 2014 - * Joerg Frohne, Texas A&M University and + * 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 @@ -109,8 +109,8 @@ namespace ElastoPlastic const double height, Triangulation<3,3> &result) { - // Assert (input.n_levels() == 1, - // ExcMessage ("The input triangulations must be coarse meshes.")); + // 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.")); @@ -124,13 +124,13 @@ namespace ElastoPlastic for (unsigned int i=0; i &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); - } + if (input.get_used_vertices()[i]) + { + const Point<2> &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); + } } } @@ -202,117 +202,117 @@ namespace ElastoPlastic { - 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)); - } + 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)); + } } @@ -337,9 +337,9 @@ namespace ElastoPlastic template struct PointHistory { - SymmetricTensor<2,dim> old_stress; - SymmetricTensor<2,dim> old_strain; - Point point; + SymmetricTensor<2,dim> old_stress; + SymmetricTensor<2,dim> old_strain; + Point point; }; @@ -381,7 +381,7 @@ namespace ElastoPlastic bool get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, - const std::vector > &point_hessian, + const std::vector > &point_hessian, Tensor<5, dim> &stress_strain_tensor_grad) const; void @@ -490,7 +490,7 @@ namespace ElastoPlastic bool ConstitutiveLaw:: get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, - const std::vector > &point_hessian, + const std::vector > &point_hessian, Tensor<5, dim> &stress_strain_tensor_grad) const { SymmetricTensor<2, dim> stress_tensor; @@ -502,53 +502,54 @@ namespace ElastoPlastic 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; - } + { + 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); } @@ -619,18 +620,18 @@ namespace ElastoPlastic 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); + // 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); } @@ -638,65 +639,65 @@ namespace ElastoPlastic 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); + // 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); } @@ -709,129 +710,129 @@ namespace ElastoPlastic 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]); - } + /* + 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} @@ -864,30 +865,30 @@ namespace ElastoPlastic class IncrementalBoundaryValues : public Function { public: - IncrementalBoundaryValues (const double present_time, - const double present_timestep); + IncrementalBoundaryValues (const double present_time, + const double present_timestep); - virtual - void - vector_value (const Point &p, - Vector &values) const; + virtual + void + vector_value (const Point &p, + Vector &values) const; - virtual - void - vector_value_list (const std::vector > &points, - std::vector > &value_list) 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; + const double velocity; + const double present_time; + const double present_timestep; }; template IncrementalBoundaryValues:: IncrementalBoundaryValues (const double present_time, - const double present_timestep) + const double present_timestep) : Function (dim), velocity (.1), @@ -900,13 +901,13 @@ namespace ElastoPlastic void IncrementalBoundaryValues:: vector_value (const Point &p, - Vector &values) const + Vector &values) const { - Assert (values.size() == dim, - ExcDimensionMismatch (values.size(), dim)); + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); - values = 0; - values(2) = -present_timestep * velocity; + values = 0; + values(2) = -present_timestep * velocity; } @@ -915,829 +916,829 @@ namespace ElastoPlastic void IncrementalBoundaryValues:: vector_value_list (const std::vector > &points, - std::vector > &value_list) const + std::vector > &value_list) const { - const unsigned int n_points = points.size(); + const unsigned int n_points = points.size(); - Assert (value_list.size() == n_points, - ExcDimensionMismatch (value_list.size(), n_points)); + 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]); + 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]); - } - - // ------------------------------------------------------------------------------- - } - + // ----------------------------- TimoshenkoBeam --------------------------------------- + /* + template + class IncrementalBoundaryForce : public Function + { + public: + IncrementalBoundaryForce (const double present_time, + const double end_time); - namespace DualFunctional - { + virtual + void vector_value (const Point &p, + Vector &values) const; - 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; - }; + 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 - class PointValuesEvaluation : public DualFunctionalBase + void + IncrementalBoundaryForce::vector_value (const Point &p, + Vector &values) const { - public: - PointValuesEvaluation (const Point &evaluation_point); + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + AssertThrow (dim == 2, ExcNotImplemented()); - virtual - void - assemble_rhs (const DoFHandler &dof_handler, - const Vector &solution, - const ConstitutiveLaw &constitutive_law, - const DoFHandler &dof_handler_dual, - Vector &rhs_dual) const; + // compute traction on the right face of Timoshenko beam problem, t_bar + double inertia_moment = (thickness*std::pow(depth,3)) / 12; - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); + double x = p(0); + double y = p(1); - protected: - const Point evaluation_point; - }; + AssertThrow(std::fabs(x-length)<1e-12, ExcNotImplemented()); + values(0) = 0; + values(1) = - shear_force/(2*inertia_moment) * ( depth*depth/4-y*y ); - template - PointValuesEvaluation:: - PointValuesEvaluation (const Point &evaluation_point) - : - evaluation_point (evaluation_point) - {} + // compute the fraction of imposed force + const double frac = present_time/end_time; + values *= frac; + } 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 + IncrementalBoundaryForce:: + vector_value_list (const std::vector > &points, + std::vector > &value_list) const { - rhs_dual.reinit (dof_handler_dual.n_dofs()); - const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex; + const unsigned int n_points = points.size(); - 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; - } + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); - AssertThrow (false, ExcEvaluationPointNotFound(evaluation_point)); + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); } template - class PointXDerivativesEvaluation : public DualFunctionalBase + class BodyForce : public ZeroFunction { public: - PointXDerivativesEvaluation (const Point &evaluation_point); + BodyForce () : ZeroFunction (dim) {} + }; + + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double end_time); virtual void - assemble_rhs (const DoFHandler &dof_handler, - const Vector &solution, - const ConstitutiveLaw &constitutive_law, - const DoFHandler &dof_handler_dual, - Vector &rhs_dual) const; + vector_value (const Point &p, + Vector &values) const; - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; - protected: - const Point evaluation_point; + private: + const double present_time, + end_time, + shear_force, + Youngs_modulus, + Poissons_ratio, + length, + depth, + thickness; }; template - PointXDerivativesEvaluation:: - PointXDerivativesEvaluation (const Point &evaluation_point) - : - evaluation_point (evaluation_point) + 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 - PointXDerivativesEvaluation:: - assemble_rhs (const DoFHandler &dof_handler, - const Vector &solution, - const ConstitutiveLaw &constitutive_law, - const DoFHandler &dof_handler_dual, - Vector &rhs_dual) const + IncrementalBoundaryValues:: + vector_value (const Point &p, + Vector &values) const { - rhs_dual.reinit (dof_handler_dual.n_dofs()); - const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex; + AssertThrow (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + AssertThrow (dim == 2, ExcNotImplemented()); - 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); + // compute exact displacement of Timoshenko beam problem, u_bar + double inertia_moment = (thickness*std::pow(depth,3)) / 12; - double total_volume = 0; + double x = p(0); + double y = p(1); - 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; + double fac = shear_force / (6*Youngs_modulus*inertia_moment); - for (unsigned int q=0; qget_dof_indices (local_dof_indices); - for (unsigned int i=0; i + 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, @@ -1752,126 +1753,126 @@ namespace ElastoPlastic class MeanDisplacementFace : public DualFunctionalBase { public: - MeanDisplacementFace (const unsigned int face_id, - const std::vector comp_mask); + 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; + 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; + 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) + const std::vector comp_mask ) + : + face_id (face_id), + comp_mask (comp_mask) { - AssertThrow(comp_mask.size() == dim, - ExcDimensionMismatch (comp_mask.size(), dim) ); + 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 + 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= 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 { public: - MeanStressFace (const unsigned int face_id, - const std::vector > &comp_stress); + 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; + 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; + 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) + const std::vector > &comp_stress ) + : + face_id (face_id), + comp_stress (comp_stress) { - AssertThrow(comp_stress.size() == dim, - ExcDimensionMismatch (comp_stress.size(), dim) ); + 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 + 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()); + AssertThrow (dim >= 2, ExcNotImplemented()); - rhs_dual.reinit (dof_handler_dual.n_dofs()); + rhs_dual.reinit (dof_handler_dual.n_dofs()); - const QGauss face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1); + 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); + 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(); + 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; + std::vector > strain_tensor(n_face_q_points); + SymmetricTensor<4, dim> stress_strain_tensor; - Vector cell_rhs (dofs_per_cell_dual); + Vector cell_rhs (dofs_per_cell_dual); - std::vector local_dof_indices (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.; + // 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; + bool evaluation_face_found = false; - typename DoFHandler::active_cell_iterator - cell_dual = dof_handler_dual.begin_active(), - endc_dual = dof_handler_dual.end(), + 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); + 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; - } + 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.reinit (cell, face); + fe_face_values_dual.reinit (cell_dual, face); - fe_face_values[displacement].get_function_symmetric_gradients(solution, - strain_tensor); + 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 i=0; i + 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); - } + 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; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i { public: - MeanStressDomain (const std::string &base_mesh, - const std::vector > &comp_stress); + 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; + 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; + 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) + 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) ); + 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 + 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()); + 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 + // 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; + const double height = 200e-3, + thickness_flange = 10e-3; - rhs_dual.reinit (dof_handler_dual.n_dofs()); + rhs_dual.reinit (dof_handler_dual.n_dofs()); - const QGauss quadrature_formula(dof_handler_dual.get_fe().tensor_degree()+1); + 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); + 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(); + 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; + std::vector > strain_tensor(n_q_points); + SymmetricTensor<4, dim> stress_strain_tensor; - Vector cell_rhs (dofs_per_cell_dual); + Vector cell_rhs (dofs_per_cell_dual); - std::vector local_dof_indices (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.; + // 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; + bool evaluation_domain_found = false; - typename DoFHandler::active_cell_iterator - cell_dual = dof_handler_dual.begin_active(), - endc_dual = dof_handler_dual.end(), + 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); + const FEValuesExtractors::Vector displacement(0); - fe_values[displacement].get_function_symmetric_gradients(solution, - strain_tensor); + 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; - for (unsigned int q_point=0; q_point - stress_phi_i = stress_strain_tensor - * fe_values_dual[displacement].symmetric_gradient(i, q_point); + 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; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i { public: - MeanStrainEnergyFace (const unsigned int face_id, - const Function &lambda_function, - const Function &mu_function ); + 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; + 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; + 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) + const Function &lambda_function, + const Function &mu_function ) + : + face_id (face_id), + lambda_function (&lambda_function), + mu_function (&mu_function) {} @@ -2220,151 +2221,151 @@ namespace ElastoPlastic void MeanStrainEnergyFace:: assemble_rhs_nonlinear (const DoFHandler &primal_dof_handler, - const Vector &primal_solution, - const DoFHandler &dof_handler, - Vector &rhs) const + 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_point= 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; iget_dof_indices (local_dof_indices); - for (unsigned int i=0; i &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); + 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); + void compute_error_DWR (Vector &estimated_error_per_cell); - ~DualSolver (); + ~DualSolver (); private: void setup_system (); void compute_dirichlet_constraints (); - void assemble_matrix (); - void assemble_rhs (); - void solve (); - void output_results (); + void assemble_matrix (); + void assemble_rhs (); + void solve (); + void output_results (); - const FESystem fe; - DoFHandler dof_handler; - const Vector solution; + const FESystem fe; + DoFHandler dof_handler; + const Vector solution; - const unsigned int fe_degree; + const unsigned int fe_degree; - const unsigned int fe_degree_dual; - FESystem fe_dual; - DoFHandler dof_handler_dual; + const unsigned int fe_degree_dual; + FESystem fe_dual; + DoFHandler dof_handler_dual; const QGauss quadrature_formula; - const QGauss face_quadrature_formula; + const QGauss face_quadrature_formula; - ConstraintMatrix constraints_hanging_nodes_dual; - ConstraintMatrix constraints_dirichlet_and_hanging_nodes_dual; + 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; + SparsityPattern sparsity_pattern_dual; + SparseMatrix system_matrix_dual; + Vector system_rhs_dual; + Vector solution_dual; const ConstitutiveLaw constitutive_law; - const SmartPointer > triangulation; + const SmartPointer > triangulation; const SmartPointer > dual_functional; - unsigned int timestep_no; + unsigned int timestep_no; std::string output_dir; const std::string base_mesh; - double present_time; - double end_time; + 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) + 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 (); + dof_handler_dual.clear (); } template void DualSolver::setup_system() { - dof_handler.distribute_dofs(fe); + 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; + 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 (); + 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); + 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); +// constraints_hanging_nodes_dual.condense (sparsity_pattern_dual); + constraints_dirichlet_and_hanging_nodes_dual.condense (sparsity_pattern_dual); - sparsity_pattern_dual.compress(); + sparsity_pattern_dual.compress(); - system_matrix_dual.reinit (sparsity_pattern_dual); + system_matrix_dual.reinit (sparsity_pattern_dual); - solution_dual.reinit (dof_handler_dual.n_dofs()); - system_rhs_dual.reinit (dof_handler_dual.n_dofs()); + solution_dual.reinit (dof_handler_dual.n_dofs()); + system_rhs_dual.reinit (dof_handler_dual.n_dofs()); } @@ -2517,64 +2518,78 @@ namespace ElastoPlastic 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()); - } + { + 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(); } @@ -2586,10 +2601,10 @@ namespace ElastoPlastic FEValues fe_values(fe, quadrature_formula, update_gradients); FEValues fe_values_dual(fe_dual, quadrature_formula, - update_values | update_gradients | update_JxW_values); + 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(); + const unsigned int n_q_points = quadrature_formula.size(); FullMatrix cell_matrix (dofs_per_cell_dual, dofs_per_cell_dual); @@ -2605,10 +2620,10 @@ namespace ElastoPlastic for (; cell_dual != endc_dual; ++cell_dual, ++cell) if (cell_dual->is_locally_owned()) { - fe_values.reinit(cell); + fe_values.reinit(cell); - fe_values_dual.reinit(cell_dual); - cell_matrix = 0; + fe_values_dual.reinit(cell_dual); + cell_matrix = 0; std::vector > strain_tensor(n_q_points); fe_values[displacement].get_function_symmetric_gradients(solution, @@ -2616,31 +2631,31 @@ namespace ElastoPlastic 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); + 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 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)); + 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); + local_dof_indices, + system_matrix_dual); } @@ -2650,54 +2665,54 @@ namespace ElastoPlastic 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); + 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); + // +++ 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); +// PreconditionSSOR<> preconditioner; +// preconditioner.initialize(system_matrix_dual, 1.2); // -// cg.solve (system_matrix_dual, solution_dual, system_rhs_dual, -// preconditioner); +// 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); + // ++++ 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); +// PreconditionJacobi<> preconditioner; +// preconditioner.initialize(system_matrix_dual, 1.0); // -// bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual, -// preconditioner); +// bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual, +// preconditioner); - // +++++++++++++++++++++++++++++++++++++++++++++++++ + // +++++++++++++++++++++++++++++++++++++++++++++++++ - constraints_dirichlet_and_hanging_nodes_dual.distribute (solution_dual); + 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"); + 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); @@ -2727,17 +2742,17 @@ namespace ElastoPlastic 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())); + 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 (); + // solve the dual problem + setup_system (); + assemble_matrix (); + assemble_rhs (); + solve (); + output_results (); - // compuate the dual weights + // compuate the dual weights Vector primal_solution (dof_handler_dual.n_dofs()); FETools::interpolate (dof_handler, solution, @@ -2747,7 +2762,7 @@ namespace ElastoPlastic ConstraintMatrix constraints_hanging_nodes; DoFTools::make_hanging_node_constraints (dof_handler, - constraints_hanging_nodes); + constraints_hanging_nodes); constraints_hanging_nodes.close(); Vector dual_weights (dof_handler_dual.n_dofs()); FETools::interpolation_difference (dof_handler_dual, @@ -2758,399 +2773,400 @@ namespace ElastoPlastic 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(); - - } + 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(); + + } } @@ -3197,7 +3213,7 @@ namespace ElastoPlastic void setup_system (); void compute_dirichlet_constraints (); void assemble_newton_system (const TrilinosWrappers::MPI::Vector &linearization_point, - const TrilinosWrappers::MPI::Vector &delta_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 (); @@ -3278,7 +3294,7 @@ namespace ElastoPlastic // 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; + 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 @@ -3354,28 +3370,28 @@ namespace ElastoPlastic struct ErrorEstimationStrategy { - enum value - { - kelly_error, - residual_error, - weighted_residual_error, - weighted_kelly_error - }; + enum value + { + kelly_error, + residual_error, + weighted_residual_error, + weighted_kelly_error + }; }; - typename ErrorEstimationStrategy::value error_estimation_strategy; + typename ErrorEstimationStrategy::value error_estimation_strategy; - Vector estimated_error_per_cell; + Vector estimated_error_per_cell; const bool transfer_solution; std::string output_dir; - TableHandler table_results, - table_results_2, - table_results_3; + TableHandler table_results, + table_results_2, + table_results_3; unsigned int current_refinement_cycle; - const double max_relative_error; - float relative_error; + const double max_relative_error; + float relative_error; const bool show_stresses; }; @@ -3406,14 +3422,14 @@ namespace ElastoPlastic " 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"); + 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."); + 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."); @@ -3428,20 +3444,20 @@ namespace ElastoPlastic "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"), + "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)"); + 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"); + 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)"); + 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"); + 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."); @@ -3504,13 +3520,13 @@ namespace ElastoPlastic strat = prm.get("error estimation strategy"); if (strat == "kelly_error") - error_estimation_strategy = ErrorEstimationStrategy::kelly_error; + error_estimation_strategy = ErrorEstimationStrategy::kelly_error; else if (strat == "residual_error") - error_estimation_strategy = ErrorEstimationStrategy::residual_error; + error_estimation_strategy = ErrorEstimationStrategy::residual_error; else if (strat == "weighted_residual_error") - error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error; + error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error; else - AssertThrow(false, ExcNotImplemented()); + AssertThrow(false, ExcNotImplemented()); output_dir = prm.get("output directory"); if (output_dir != "" && *(output_dir.rbegin()) != '/') @@ -3550,468 +3566,485 @@ namespace ElastoPlastic 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()); - } + { + 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; + << triangulation.n_active_cells() + << std::endl; } @@ -4032,13 +4065,13 @@ namespace ElastoPlastic 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"); { 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; + << dof_handler.n_dofs() + << std::endl; locally_owned_dofs = dof_handler.locally_owned_dofs(); locally_relevant_dofs.clear(); @@ -4066,9 +4099,9 @@ namespace ElastoPlastic { TimerOutput::Scope t(computing_timer, "Setup: vectors"); if (timestep_no==1 || current_refinement_cycle!=0) - { - solution.reinit(locally_relevant_dofs, mpi_communicator); - } + { + 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); @@ -4119,64 +4152,78 @@ namespace ElastoPlastic 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()); - } + { + 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(); @@ -4197,21 +4244,23 @@ namespace ElastoPlastic void ElastoPlasticProblem:: assemble_newton_system (const TrilinosWrappers::MPI::Vector &linearization_point, - const TrilinosWrappers::MPI::Vector &delta_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; - } + { + 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 | @@ -4225,14 +4274,14 @@ namespace ElastoPlastic 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::BodyForce body_force; + std::vector > body_force_values(n_q_points, + Vector(dim)); const EquationData:: - IncrementalBoundaryForce boundary_force(present_time, end_time); + IncrementalBoundaryForce boundary_force(present_time, end_time); std::vector > boundary_force_values(n_face_q_points, - Vector(dim)); + Vector(dim)); FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); Vector cell_rhs(dofs_per_cell); @@ -4256,7 +4305,7 @@ namespace ElastoPlastic cell_rhs = 0; fe_values[displacement].get_function_symmetric_gradients(delta_linearization_point, - incremental_strain_tensor); + incremental_strain_tensor); // For assembling the local right hand side contributions, we need // to access the prior linearized stress value in this quadrature @@ -4266,13 +4315,13 @@ namespace ElastoPlastic // 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()); + = 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: @@ -4281,9 +4330,9 @@ namespace ElastoPlastic 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<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; @@ -4293,9 +4342,9 @@ namespace ElastoPlastic 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]; - } + { + rhs_values_body_force[i] = body_force_values[q_point][i]; + } for (unsigned int i = 0; i < dofs_per_cell; ++i) { @@ -4326,11 +4375,11 @@ namespace ElastoPlastic * fe_values.JxW(q_point)); cell_rhs(i) += ( - ( stress_phi_i - * incremental_strain_tensor[q_point] ) + ( stress_phi_i + * incremental_strain_tensor[q_point] ) - ( ( stress_strain_tensor - * fe_values[displacement].symmetric_gradient(i, q_point)) + * fe_values[displacement].symmetric_gradient(i, q_point)) * tmp_strain_tensor_qpoint ) + ( fe_values[displacement].value(i, q_point) @@ -4354,9 +4403,9 @@ namespace ElastoPlastic { Tensor<1, dim> rhs_values; for (unsigned int i = 0; i < dim; ++i) - { - rhs_values[i] = boundary_force_values[q_point][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 @@ -4366,10 +4415,10 @@ namespace ElastoPlastic 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); + local_dof_indices, + newton_matrix, + newton_rhs, + true); } @@ -4410,15 +4459,17 @@ namespace ElastoPlastic { 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; - } + { + 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 | @@ -4432,14 +4483,14 @@ namespace ElastoPlastic 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::BodyForce body_force; + std::vector > body_force_values(n_q_points, + Vector(dim)); const EquationData:: - IncrementalBoundaryForce boundary_force(present_time, end_time); + IncrementalBoundaryForce boundary_force(present_time, end_time); std::vector > boundary_force_values(n_face_q_points, - Vector(dim)); + Vector(dim)); Vector cell_rhs(dofs_per_cell); @@ -4466,7 +4517,7 @@ namespace ElastoPlastic strain_tensors); body_force.vector_value_list(fe_values.get_quadrature_points(), - body_force_values); + body_force_values); for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { @@ -4479,20 +4530,20 @@ namespace ElastoPlastic 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]; - } + { + 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); + * 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; @@ -4516,9 +4567,9 @@ namespace ElastoPlastic { Tensor<1, dim> rhs_values; for (unsigned int i = 0; i < dim; ++i) - { - rhs_values[i] = boundary_force_values[q_point][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)); @@ -4527,8 +4578,8 @@ namespace ElastoPlastic cell->get_dof_indices(local_dof_indices); constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_rhs, - local_dof_indices, - newton_rhs_residual); + local_dof_indices, + newton_rhs_residual); } @@ -4688,13 +4739,13 @@ namespace ElastoPlastic 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); + 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; + previous_disp_norm = 0; const double correct_sigma = sigma_0; @@ -4715,14 +4766,14 @@ namespace ElastoPlastic pcout << " Newton iteration " << newton_step << std::endl; pcout << " Assembling system... " << std::endl; - newton_matrix = 0; - newton_rhs = 0; + newton_matrix = 0; + newton_rhs = 0; newton_rhs_residual = 0; tmp_solution = solution; tmp_solution += incremental_displacement; assemble_newton_system(tmp_solution, - incremental_displacement); + incremental_displacement); pcout << " Solving system... " << std::endl; solve_newton_system(); @@ -4760,8 +4811,8 @@ namespace ElastoPlastic || (!transfer_solution && newton_step == 2)) { - tmp_solution = solution; - tmp_solution += incremental_displacement; + tmp_solution = solution; + tmp_solution += incremental_displacement; compute_nonlinear_residual(tmp_solution); old_solution = incremental_displacement; @@ -4817,16 +4868,16 @@ namespace ElastoPlastic // the iteration on the current mesh: // if (residual_norm < 1e-10) if (residual_norm < 1e-7) - break; + break; pcout << " difference of two consecutive incremental displacement l2 norm : " - << std::abs(disp_norm - previous_disp_norm) << std::endl; + << 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; - } + (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; @@ -4840,156 +4891,163 @@ namespace ElastoPlastic 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; - } - } + 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); - DualFunctional::MeanStressFace dual_functional(face_id, comp_stress); + } + else if (error_estimation_strategy == ErrorEstimationStrategy::weighted_residual_error) + { + // make a non-parallel copy of tmp_solution + Vector copy_solution(tmp_solution); - DualSolver dual_solver(triangulation, fe, - copy_solution, - constitutive_law, dual_functional, - timestep_no, output_dir, base_mesh, - present_time, end_time); + // the dual function definition (it should be defined previously, e.g. input file) + if (base_mesh == "Timoshenko beam") + { + double length = .48, + depth = .12; - dual_solver.compute_error_DWR (estimated_error_per_cell); + const Point evaluation_point(length, -depth/2); - }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::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); + 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); - 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); - 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; - }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; + DualFunctional::MeanStressFace dual_functional(face_id, comp_stress); - const Point evaluation_point(0, height/2, length); + // ......................................... - 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); - // ----------------------------------------------------------- - // 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 + dual_solver.compute_error_DWR (estimated_error_per_cell); - 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); + } + 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); + 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); + dual_solver.compute_error_DWR (estimated_error_per_cell); - }else - { - AssertThrow(false, ExcNotImplemented()); - } + } + else + { + AssertThrow(false, ExcNotImplemented()); + } - }else - { - AssertThrow(false, ExcNotImplemented()); - } + } + else + { + AssertThrow(false, ExcNotImplemented()); + } - relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm(); + relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm(); - pcout << "Estimated relative error = " << relative_error << std::endl; + pcout << "Estimated relative error = " << relative_error << std::endl; } @@ -4997,382 +5055,383 @@ namespace ElastoPlastic 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(); - - } + 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(); + + } } @@ -5390,81 +5449,81 @@ namespace ElastoPlastic 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 + // --------------------------------------------------------------- + // 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 @@ -5475,10 +5534,10 @@ namespace ElastoPlastic } else { - const double refine_fraction_cells = .3, - coarsen_fraction_cells = .03; -// const double refine_fraction_cells = .1, -// coarsen_fraction_cells = .3; + 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, @@ -5489,40 +5548,40 @@ namespace ElastoPlastic triangulation.prepare_coarsening_and_refinement(); parallel::distributed::SolutionTransfer solution_transfer(dof_handler); + TrilinosWrappers::MPI::Vector> solution_transfer(dof_handler); solution_transfer.prepare_for_coarsening_and_refinement(solution); parallel::distributed::SolutionTransfer incremental_displacement_transfer(dof_handler); + TrilinosWrappers::MPI::Vector> incremental_displacement_transfer(dof_handler); if (transfer_solution) - incremental_displacement_transfer.prepare_for_coarsening_and_refinement(incremental_displacement); + 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_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]); - } + { + 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]); - } + { + 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_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]); - } + { + 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]); - } + { + history_strain_field_transfer2.prepare_for_coarsening_and_refinement(history_strain_field[2]); + } triangulation.execute_coarsening_and_refinement(); pcout << " Number of active cells: " @@ -5553,20 +5612,20 @@ namespace ElastoPlastic 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 > 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_transfer2.interpolate(history_stress_field[2], distributed_history_stress_field[2]); + } history_stress_field = distributed_history_stress_field; @@ -5574,20 +5633,20 @@ namespace ElastoPlastic 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]); - } + { + 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_transfer2.interpolate(history_strain_field[2], distributed_history_strain_field[2]); + } history_strain_field = distributed_history_strain_field; @@ -5598,49 +5657,49 @@ namespace ElastoPlastic // new mesh. The following code will do that: FullMatrix dof_to_qpoint_matrix (quadrature_formula.size(), - history_fe.dofs_per_cell); + history_fe.dofs_per_cell); FETools::compute_interpolation_to_quadrature_points_matrix - (history_fe, - quadrature_formula, - dof_to_qpoint_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; qis_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; q 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()); + // 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} @@ -5776,69 +5835,69 @@ namespace ElastoPlastic 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; qFEValues 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; q::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()); - } + 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); @@ -5953,1063 +6012,1066 @@ namespace ElastoPlastic 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.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)); - // - data_out.build_patches(); + break; + } - // 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)); + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; - std::ofstream output_vtu((filename + ".vtu").c_str()); - data_out.write_vtu(output_vtu); - pcout << output_dir + filename_base << ".pvtu" << 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); - if (this_mpi_process == 0) + 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") { - 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"); + const double pressure (0.6*2.4e8), + inner_radius (.1); +// const double pressure (1.94e8), +// inner_radius (.1); - 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); + // 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))); - // produce eps files for mesh illustration - std::ofstream output_eps((filename + ".eps").c_str()); - GridOut grid_out; - grid_out.write_eps(triangulation, output_eps); - } + const Point point_A(inner_radius, 0.); + Vector disp_A(dim); - // Extrapolate the stresses from Gauss point to the nodes - SymmetricTensor<2, dim> stress_at_qpoint; + // make a non-parallel copy of solution + Vector copy_solution(solution); - 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:: + PointValuesEvaluation point_values_evaluation(point_A); - Vector 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); + point_values_evaluation.compute (dof_handler, copy_solution, disp_A); - 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); + 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)); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(), - dg_cell = history_dof_handler.begin_active(); + // 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()); - const FEValuesExtractors::Vector displacement(0); + // 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 - 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); - - - } + FEValues fe_values (fe, quadrature_formula, update_quadrature_points | + update_values | update_gradients); - // 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); + const unsigned int n_q_points = quadrature_formula.size(); - AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(), - ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices())); + std::vector > strain_tensor(n_q_points); + SymmetricTensor<4, dim> stress_strain_tensor; + Tensor<2, dim> stress_at_qpoint; - std::vector< std::vector< Vector > > - history_stress_on_vertices (dim, std::vector< Vector >(dim)); - for (unsigned int i=0; i 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_on_vertices (dof_handler_1.n_dofs()), - counter_on_vertices (dof_handler_1.n_dofs()); - VM_stress_on_vertices = 0; - counter_on_vertices = 0; + 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); - 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::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + dg_cell = history_dof_handler.begin_active(); - for (unsigned int id=0; idis_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 data_out; - data_out.attach_dof_handler (history_dof_handler); + // transform the stress from the Cartesian coordinate to the polar coordinate + const Point 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); - 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"); + // stress in polar coordinate + stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix); - 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"); - } + for (unsigned int i=0; iset_dof_values (local_history_stress_fe_values[i][j], + history_stress_field[i][j]); + } - 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"); + { + DataOut data_out; + data_out.attach_dof_handler (history_dof_handler); - 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); - } + 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); - { - DataOut data_out; - data_out.attach_dof_handler (dof_handler_1); + 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; - 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 (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"); - 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"); - } + std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str()); + data_out.write_pvtu_record(pvtu_master_output, filenames); - data_out.build_patches (); + std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str()); + data_out.write_visit_record(visit_master_output, filenames); + } - 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; + // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + // 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 filenames; - for (unsigned int i = 0; i < n_mpi_processes; ++i) - filenames.push_back(filename_base_stress + "-" + - Utilities::int_to_string(i, 4) + - ".vtu"); + Vector counter_on_vertices (dof_handler_1.n_dofs()); + counter_on_vertices = 0; - std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str()); - data_out.write_pvtu_record(pvtu_master_output, filenames); + 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()) + { - std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str()); - data_out.write_visit_record(visit_master_output, filenames); - } + 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; - magnified_solution *= -1; - move_mesh(magnified_solution); + for (unsigned int i=0; i 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)); + for (unsigned int id=0; id data_out; + data_out.attach_dof_handler (dof_handler_1); - 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); + 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"); - std::cout << "Exact displacement at the point (" << intersted_point[0] - << ", " << intersted_point[1] << ") is " - << "(" << vertex_exact_displacement[0] - << ", " << vertex_exact_displacement[1] << ").\n\n"; + data_out.build_patches (); - }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); + 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)); - // 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; + 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"); - 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; - } + 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); + } - 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); + if ( std::abs( (present_time/end_time)*(pressure/sigma_0) - 0.6 ) < + .501*(present_timestep/end_time)*(pressure/sigma_0) ) + { - vertices_found[dof_1_vertex] = true; - } - } + // 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; - AssertThrow(evaluation_face_found, ExcInternalError()); + 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; + } + } + + } + } - // 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; + AssertThrow(evaluation_face_found, ExcInternalError()); - std::vector vertices_found_2 (dof_handler_1.n_dofs(), false); - unsigned int no_vertices_found = 0; + // table_results_3: presenting the mean stress_rr of the nodes on the inner radius + const unsigned int face_id_2 = 0; - evaluation_face_found = false; + Tensor<2, dim> stress_node, + mean_stress_polar; + mean_stress_polar = 0; - 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; - } + std::vector vertices_found_2 (dof_handler_1.n_dofs(), false); + unsigned int no_vertices_found = 0; + evaluation_face_found = false; - 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; iris_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); + // 3d and plane strain + const double mu (e_modulus / (2 * (1 + nu))); - // make a non-parallel copy of solution - Vector copy_solution(solution); + // table_results: Demonstrates the result of displacement at the top left corner versus imposed tension + /* + { + const Point point_C(0., height); + Vector disp_C(dim); - typename Evaluation::PointValuesEvaluation:: - PointValuesEvaluation point_values_evaluation(point_C); + // make a non-parallel copy of solution + Vector copy_solution(solution); - point_values_evaluation.compute (dof_handler, copy_solution, disp_C); + typename Evaluation::PointValuesEvaluation:: + PointValuesEvaluation point_values_evaluation(point_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)); - } - */ + point_values_evaluation.compute (dof_handler, copy_solution, disp_C); - // compute average sigma_yy on the bottom edge - double stress_yy_av; - { - stress_yy_av = 0; - const unsigned int face_id = 1; + 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)); + } + */ - std::vector vertices_found (dof_handler_1.n_dofs(), false); - unsigned int no_vertices_in_face = 0; + // compute average sigma_yy on the bottom edge + double stress_yy_av; + { + stress_yy_av = 0; + const unsigned int face_id = 1; - bool evaluation_face_found = false; + std::vector vertices_found (dof_handler_1.n_dofs(), false); + unsigned int no_vertices_in_face = 0; - 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; - } + 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; + } + } + + } + } - 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; - } + } - AssertThrow(evaluation_face_found, ExcInternalError()); + // table_results_2: Demonstrate the stress_yy on the nodes of bottom edge - stress_yy_av /= no_vertices_in_face; +// 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); - // table_results_2: Demonstrate the stress_yy on the nodes of bottom edge + bool evaluation_face_found = false; -// 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()); + 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; + } + } + + } + } - } + } - // 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; + AssertThrow(evaluation_face_found, ExcInternalError()); - // 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); + // 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; - SymmetricTensor<2, dim> strain_at_qpoint; + // 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); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(), - dg_cell = history_dof_handler.begin_active(); + Vector local_strain_yy_values_at_qpoints (quadrature_formula.size()), + local_strain_yy_fe_values (history_fe.dofs_per_cell); - bool cell_found = false; + SymmetricTensor<2, dim> strain_at_qpoint; - 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()); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + dg_cell = history_dof_handler.begin_active(); - // Then loop over the quadrature points of this cell: - for (unsigned int q=0; qis_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); - }else if (base_mesh == "Cantiliver_beam_3d") - { - const double pressure (6e6), - length (.7), - height (200e-3); + // make a non-parallel copy of solution + Vector copy_solution(solution); - // table_results: Demonstrates the result of displacement at the top front point, Point A - { - // Quantity of interest: - // displacement at Point A (x=0, y=height/2, z=length) + typename Evaluation::PointValuesEvaluation:: + PointValuesEvaluation point_values_evaluation(point_A); - const Point point_A(0, height/2, length); - Vector disp_A(dim); + point_values_evaluation.compute (dof_handler, copy_solution, disp_A); - // make a non-parallel copy of solution - Vector copy_solution(solution); + 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)); + } - typename Evaluation::PointValuesEvaluation:: - PointValuesEvaluation point_values_evaluation(point_A); + { + // 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 - point_values_evaluation.compute (dof_handler, copy_solution, disp_A); + double VM_stress_max (0); + Point point_max; - 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)); - } + SymmetricTensor<2, dim> stress_at_qpoint; - { - // 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 + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); - double VM_stress_max (0); - Point point_max; + const FEValuesExtractors::Vector displacement(0); - SymmetricTensor<2, dim> stress_at_qpoint; + 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::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + const double VM_stress = Evaluation::get_von_Mises_stress(stress_at_qpoint); + if (VM_stress > VM_stress_max) + { + VM_stress_max = VM_stress; + point_max = local_quadrature_points_history[q].point; + } - 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]); + 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]); - } + } - } + } } @@ -7029,148 +7091,149 @@ namespace ElastoPlastic 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); - } + 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); + } } } -- 2.39.5