From d31d25c443a9b7eaee60b80b8cd27c7a1f218856 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 23 Sep 2013 21:15:23 +0000 Subject: [PATCH] More small rewrites. git-svn-id: https://svn.dealii.org/trunk@30901 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 456 +++++++++++++--------------- 1 file changed, 203 insertions(+), 253 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 275c8c3623..5c6f6a1bcb 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -94,7 +94,7 @@ namespace Step42 // // The data which we read in by the // function read_obstacle () from the file -// "obstacle_file.pbm" will be stored +// will be stored // in a double std::vector named // obstacle_data. // This vector composes the base @@ -103,142 +103,57 @@ namespace Step42 // This will be done by obstacle_function (). // // In the function run () of the class -// PlasticityContactProblem we create -// an object of the class Input which will +// PlasticityContactProblem we create +// an object of the this class which will // be used in the class Obstacle to // supply the obstacle function in // update_solution_and_constraints () of // the class PlasticityContactProblem. - +// +// The hx,hy variables denote the spacing between +// pixels in $x$ and $y$ directions. nx,ny are +// the numbers of pixels in each of these directions. +// get_value() returns the value of the +// image at a given location, interpolated from the adjacent +// pixel values. template class Input { public: - Input (const std::string &name) - : - pcout(std::cout, - (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)), - obstacle_data(0), - hx(0), - hy(0), - nx(0), - ny(0) - { - read_obstacle(name); - } - - double hv (const int i, - const int j); + Input (const std::string &name); double - obstacle_function (const double x, - const double y); - - void - read_obstacle (const std::string name); + get_value (const double x, + const double y); private: - ConditionalOStream pcout; std::vector obstacle_data; double hx, hy; int nx, ny; - }; - -// This function is used in obstacle_function () -// to provide the proper value of the obstacle. - template - double - Input::hv (const int i, - const int j) - { - assert(i >= 0 && i < nx); - assert(j >= 0 && j < ny); - return obstacle_data[nx * (ny - 1 - j) + i]; // i indiziert x-werte, j indiziert y-werte - } -// obstacle_function () calculates the bilinear interpolated -// value in the point (x,y). - template - double - Input::obstacle_function (const double x, - const double y) - { - int ix = (int) (x / hx); - int iy = (int) (y / hy); - - if (ix < 0) - ix = 0; - - if (iy < 0) - iy = 0; - - if (ix >= nx - 1) - ix = nx - 2; - - if (iy >= ny - 1) - iy = ny - 2; - - double val = 0.0; - { - FullMatrix H(4, 4); - Vector X(4); - Vector b(4); - - double xx, yy; - - xx = ix * hx; - yy = iy * hy; - H(0, 0) = xx; - H(0, 1) = yy; - H(0, 2) = xx * yy; - H(0, 3) = 1.0; - b(0) = hv(ix, iy); - - xx = (ix + 1) * hx; - yy = iy * hy; - H(1, 0) = xx; - H(1, 1) = yy; - H(1, 2) = xx * yy; - H(1, 3) = 1.0; - b(1) = hv(ix + 1, iy); - - xx = (ix + 1) * hx; - yy = (iy + 1) * hy; - H(2, 0) = xx; - H(2, 1) = yy; - H(2, 2) = xx * yy; - H(2, 3) = 1.0; - b(2) = hv(ix + 1, iy + 1); - - xx = ix * hx; - yy = (iy + 1) * hy; - H(3, 0) = xx; - H(3, 1) = yy; - H(3, 2) = xx * yy; - H(3, 3) = 1.0; - b(3) = hv(ix, iy + 1); - - H.gauss_jordan(); - H.vmult(X, b); - - val = X(0) * x + X(1) * y + X(2) * x * y + X(3); - } + double get_pixel_value (const int i, + const int j); + }; - return val; - } -// As mentioned above this function reads in the -// obstacle data and stores them in the std::vector -// obstacle_data. It will be used only in run (). - template - void - Input::read_obstacle (const std::string name) + // The constructor of this class reads in the data that describes + // the obstacle from the given file name. + template + Input::Input(const std::string &name) + : + obstacle_data(0), + hx(0), + hy(0), + nx(0), + ny(0) { std::ifstream f(name.c_str()); std::string temp; f >> temp >> nx >> ny; - assert(nx > 0 && ny > 0); + + AssertThrow(nx > 0 && ny > 0, + ExcMessage ("Invalid file format.")); for (int k = 0; k < nx * ny; k++) { @@ -250,14 +165,80 @@ namespace Step42 hx = 1.0 / (nx - 1); hy = 1.0 / (ny - 1); - pcout << "Resolution of the scanned obstacle picture: " << nx << " x " - << ny << std::endl; + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + std::cout << "Resolution of the scanned obstacle picture: " << nx << " x " << ny + << std::endl; + } + + template + double + Input::get_pixel_value (const int i, + const int j) + { + assert(i >= 0 && i < nx); + assert(j >= 0 && j < ny); + return obstacle_data[nx * (ny - 1 - j) + i]; + } + + + template + double + Input::get_value (const double x, + const double y) + { + const int ix = std::min(std::max((int) (x / hx), 0), nx-2); + const int iy = std::min(std::max((int) (y / hy), 0), ny-2); + + FullMatrix H(4, 4); + Vector X(4); + Vector b(4); + + double xx, yy; + + xx = ix * hx; + yy = iy * hy; + H(0, 0) = xx; + H(0, 1) = yy; + H(0, 2) = xx * yy; + H(0, 3) = 1.0; + b(0) = get_value(ix, iy); + + xx = (ix + 1) * hx; + yy = iy * hy; + H(1, 0) = xx; + H(1, 1) = yy; + H(1, 2) = xx * yy; + H(1, 3) = 1.0; + b(1) = get_pixel_value(ix + 1, iy); + + xx = (ix + 1) * hx; + yy = (iy + 1) * hy; + H(2, 0) = xx; + H(2, 1) = yy; + H(2, 2) = xx * yy; + H(2, 3) = 1.0; + b(2) = get_pixel_value(ix + 1, iy + 1); + + xx = ix * hx; + yy = (iy + 1) * hy; + H(3, 0) = xx; + H(3, 1) = yy; + H(3, 2) = xx * yy; + H(3, 3) = 1.0; + b(3) = get_pixel_value(ix, iy + 1); + + H.gauss_jordan(); + H.vmult(X, b); + + return (X(0) * x + X(1) * y + X(2) * x * y + X(3)); } + // @sect3{The ConstitutiveLaw class template} // This class provides an interface -// for a constitutive law. In this +// for a constitutive law, i.e., for the relationship between strain +// $\varepsilon(\mathbf u)$ and stress $\sigma$. In this // example we are using an elastoplastic // material behavior with linear, // isotropic hardening. @@ -273,18 +254,13 @@ namespace Step42 const double _gamma); bool - plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor) const; + get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor) const; void - linearized_plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor_linearized, - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor) const; - - SymmetricTensor<2, dim> - get_strain (const FEValues &fe_values, - const unsigned int shape_func, - const unsigned int q_point) const; + get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor_linearized, + SymmetricTensor<4, dim> &stress_strain_tensor) const; void set_sigma_0 (double sigma_zero) @@ -293,64 +269,47 @@ namespace Step42 } private: - const double E; - const double nu; double sigma_0; const double gamma; - const double mu; const double kappa; + const double mu; const SymmetricTensor<4, dim> stress_strain_tensor_kappa; const SymmetricTensor<4, dim> stress_strain_tensor_mu; }; // The constructor of the ConstitutiveLaw class sets the -// required material parameter for our deformable body: -// E -> elastic modulus -// nu -> Passion's number -// sigma_0 -> yield stress -// gamma -> hardening parameter. -// Also it supplies the stress strain tensor of forth order -// of the volumetric and deviator part. For further details -// see the documentation above. +// required material parameter for our deformable body. Material +// parameters for elastic isotropic media can be defined in a +// variety of ways, such as the pair $E, \nu$ (elastic modulus and +// Poisson's number), using the Lame parameters $\lambda,mu$ or +// several other commonly used conventions. Here, the constructor takes a description of material parameters in the form of $E,\nu$, but since this turns out to these are not the coefficients that appear in the equations of the plastic projector, we immediately convert them into the more suitable set $\kappa,\mu$ of bulk and shear moduli. +// In addition, the constructor takes $\sigma_0$ (the yield stress absent any plastic strain) and +// $\gamma$ (the hardening parameter) as arguments. In this constructor, we also compute the two principal components of the +// stress-strain relation and its linearization. template - ConstitutiveLaw::ConstitutiveLaw (double _E, - double _nu, - double _sigma_0, - double _gamma) + ConstitutiveLaw::ConstitutiveLaw (double E, + double nu, + double sigma_0, + double gamma) : - E(_E), - nu(_nu), - sigma_0(_sigma_0), - gamma(_gamma), - mu (E / (2 * (1 + nu))), + sigma_0(sigma_0), + gamma(gamma), kappa (E / (3 * (1 - 2 * nu))), + mu (E / (2 * (1 + nu))), stress_strain_tensor_kappa (kappa - * outer_product(unit_symmetric_tensor(), - unit_symmetric_tensor())), + * outer_product(unit_symmetric_tensor(), + unit_symmetric_tensor())), stress_strain_tensor_mu (2 * mu - * (identity_tensor() - - outer_product(unit_symmetric_tensor(), - unit_symmetric_tensor()) / 3.0)) + * (identity_tensor() + - outer_product(unit_symmetric_tensor(), + unit_symmetric_tensor()) / 3.0)) {} -// @sect4{ConstitutiveLaw::ConstitutiveLaw} -// Calculates the strain $\varepsilon(\varphi)=\dfrac{1}{2}\left(\nabla\varphi + \nabla\varphi^T$ -// for the shape functions $\varphi$. - template - SymmetricTensor<2, dim> - ConstitutiveLaw::get_strain (const FEValues &fe_values, - const unsigned int shape_func, - const unsigned int q_point) const - { - const FEValuesExtractors::Vector displacement(0); - return fe_values[displacement].symmetric_gradient(shape_func, q_point); - } +// @sect4{ConstitutiveLaw::get_stress_strain_tensor} -// @sect4{ConstitutiveLaw::plast_linear_hardening} - -// This is the implemented constitutive law. It projects the +// This is the principal component of the constitutive law. It projects the // deviatoric part of the stresses in a quadrature point back to // the yield stress (i.e., the original yield stress $\sigma_0$ plus // the term that describes linear isotropic hardening). @@ -364,8 +323,8 @@ namespace Step42 template bool ConstitutiveLaw:: - plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor) const + get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor) const { Assert (dim == 3, ExcNotImplemented()); @@ -388,24 +347,25 @@ namespace Step42 return (deviator_stress_tensor_norm > sigma_0); } -// @sect4{ConstitutiveLaw::linearized_plast_linear_hardening} + +// @sect4{ConstitutiveLaw::get_linearized_stress_strain_tensors} // This function returns the linearized stress strain tensor, linearized // around the solution $u^{i-1}$ of the previous Newton step $i-1$. -// The parameter strain_tensor $\varepsilon(u^{i-1})$ is calculated -// by $u^{i-1}$. It contains the derivative of the nonlinear -// constitutive law. As the result this function returns -// the stress_strain_tensor of the nonlinear problem as well as -// the stress_strain_tensor_linearized of the linearized problem. +// The parameter strain_tensor (commonly denoted $\varepsilon(u^{i-1})$) must be passed as an argument, +// and serves as the linearization point. The function returns the derivative of the nonlinear +// constitutive law in +// the variable stress_strain_tensor, as well as +// the stress-strain tensor of the linearized problem in stress_strain_tensor_linearized. // See // PlasticityContactProblem::assemble_nl_system(TrilinosWrappers::MPI::Vector &u) // where this function is used. template void ConstitutiveLaw:: - linearized_plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor_linearized, - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor) const + get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor, + SymmetricTensor<4, dim> &stress_strain_tensor_linearized, + SymmetricTensor<4, dim> &stress_strain_tensor) const { Assert (dim == 3, ExcNotImplemented()); @@ -434,98 +394,89 @@ namespace Step42 stress_strain_tensor_linearized += stress_strain_tensor_kappa; } + //

Equation data: right hand side and boundary values

+ // + // The following should be relatively standard. We need classes for + // the right hand side forcing term (which we here choose to be zero) + // and boundary values on those part of the boundary that are not part + // of the contact surface (also chosen to be zero here). namespace EquationData { -// It possible to apply an additional body force -// but in here it is set to zero. template class RightHandSide : public Function { public: - RightHandSide () - : - Function(dim) - { - } + RightHandSide (); - virtual double - value ( - const Point &p, const unsigned int component = 0) const; + virtual + double value (const Point &p, + const unsigned int component = 0) const; - virtual void - vector_value ( - const Point &p, Vector &values) const; + virtual + void vector_value (const Point &p, + Vector &values) const; }; template - double - RightHandSide::value ( - const Point &p, const unsigned int component) const - { - double return_value = 0.0; + RightHandSide::RightHandSide () + : + Function(dim) + {} - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - return_value = 0.0; - return return_value; + template + double + RightHandSide::value (const Point &, + const unsigned int) const + { + return 0.; } template void - RightHandSide::vector_value ( - const Point &p, Vector &values) const + RightHandSide::vector_value (const Point &p, + Vector &values) const { for (unsigned int c = 0; c < this->n_components; ++c) values(c) = RightHandSide::value(p, c); } -// This function class is used to describe the prescribed displacements -// at the boundary. But again we set this to zero. + + template class BoundaryValues : public Function { public: - BoundaryValues () - : - Function(dim) - { - } - ; + BoundaryValues (); - virtual double - value ( - const Point &p, const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - virtual void - vector_value ( - const Point &p, Vector &values) const; + virtual + void vector_value (const Point &p, + Vector &values) const; }; + template - double - BoundaryValues::value ( - const Point &p, const unsigned int component) const - { - double return_value = 0; + BoundaryValues::BoundaryValues () + : + Function(dim) + {} - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - return_value = 0.0; - return return_value; + template + double + BoundaryValues::value (const Point &, + const unsigned int) const + { + return 0.; } template void - BoundaryValues::vector_value ( - const Point &p, Vector &values) const + BoundaryValues::vector_value (const Point &p, + Vector &values) const { for (unsigned int c = 0; c < this->n_components; ++c) values(c) = BoundaryValues::value(p, c); @@ -546,7 +497,7 @@ namespace Step42 bool _use_read_obstacle, double z_max_domain) : Function(dim), - input_obstacle_copy(_input), + input_obstacle(_input), use_read_obstacle(_use_read_obstacle), z_max_domain(z_max_domain) { @@ -561,7 +512,7 @@ namespace Step42 const Point &p, Vector &values) const; private: - std_cxx1x::shared_ptr > const &input_obstacle_copy; + const std_cxx1x::shared_ptr > &input_obstacle; bool use_read_obstacle; double z_max_domain; }; @@ -581,7 +532,7 @@ namespace Step42 { if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0) return z_max_domain + 0.999 - - input_obstacle_copy->obstacle_function(p(0), p(1)); + - input_obstacle->get_value(p(0), p(1)); else return 10000.0; } @@ -1039,19 +990,19 @@ namespace Step42 SymmetricTensor<4, dim> stress_strain_tensor; SymmetricTensor<2, dim> stress_tensor; - plast_lin_hard->linearized_plast_linear_hardening( - stress_strain_tensor_linearized, stress_strain_tensor, - strain_tensor[q_point]); + plast_lin_hard->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; ++i) { stress_tensor = stress_strain_tensor_linearized - * plast_lin_hard->get_strain(fe_values, i, q_point); + * fe_values[displacement].symmetric_gradient(i, q_point); for (unsigned int j = 0; j < dofs_per_cell; ++j) { cell_matrix(i, j) += (stress_tensor - * plast_lin_hard->get_strain(fe_values, j, q_point) + * fe_values[displacement].symmetric_gradient(i, q_point) * fe_values.JxW(q_point)); } @@ -1062,7 +1013,7 @@ namespace Step42 // the residual part a(v^i;v) of the rhs cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor - * plast_lin_hard->get_strain(fe_values, i, q_point) + * fe_values[displacement].symmetric_gradient(i, q_point) * fe_values.JxW(q_point)); // the residual part F(v) of the rhs @@ -1169,22 +1120,21 @@ namespace Step42 SymmetricTensor<2, dim> stress_tensor; const bool q_point_is_plastic - = plast_lin_hard->plast_linear_hardening(stress_strain_tensor, - strain_tensor[q_point]); + = plast_lin_hard->get_stress_strain_tensor(strain_tensor[q_point], + stress_strain_tensor); if (q_point_is_plastic) - { - ++plast_points; + { + ++plast_points; ++cell_constitution(cell_number); - } + } else - ++elast_points; + ++elast_points; for (unsigned int i = 0; i < dofs_per_cell; ++i) { cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor - * //(stress_tensor) * - plast_lin_hard->get_strain(fe_values, i, q_point) + * fe_values[displacement].symmetric_gradient(i, q_point) * fe_values.JxW(q_point)); Tensor<1, dim> rhs_values; @@ -1203,7 +1153,7 @@ namespace Step42 fe_values_face.reinit(cell, face); right_hand_side.vector_value_list(fe_values_face.get_quadrature_points(), - right_hand_side_values_face); + right_hand_side_values_face); for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) @@ -1241,9 +1191,9 @@ namespace Step42 // constraints_hanging_nodes.condense(system_rhs_lambda); const unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, - mpi_communicator); + mpi_communicator); const unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, - mpi_communicator); + mpi_communicator); pcout << " Number of elastic quadrature points: " << sum_elast_points << " and plastic quadrature points: " << sum_plast_points << std::endl; -- 2.39.5