From 1551a83596551484d992375797b9c682af49aa3b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 23 Sep 2013 13:41:10 +0000 Subject: [PATCH] Use standard indentation scheme. git-svn-id: https://svn.dealii.org/trunk@30893 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 3068 +++++++++++++-------------- 1 file changed, 1534 insertions(+), 1534 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 37f3257cf9..f1fd4d9db3 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -111,152 +111,152 @@ namespace Step42 // the class PlasticityContactProblem. template - class Input + class Input + { + public: + Input (const std::string &name) + : + mpi_communicator(MPI_COMM_WORLD), + pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), + obstacle_data(0), + hx(0), + hy(0), + nx(0), + ny(0) { - public: - Input (const std::string &name) - : - mpi_communicator(MPI_COMM_WORLD), - pcout(std::cout, - (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), - obstacle_data(0), - hx(0), - hy(0), - nx(0), - ny(0) - { - read_obstacle(name); - } + read_obstacle(name); + } - double - hv ( - int i, int j); + double + hv ( + int i, int j); - double - obstacle_function (const double x, - const double y); + double + obstacle_function (const double x, + const double y); - void - read_obstacle (const std::string name); + void + read_obstacle (const std::string name); - private: - MPI_Comm mpi_communicator; - ConditionalOStream pcout; - std::vector obstacle_data; - double hx, hy; - int nx, ny; - }; + private: + MPI_Comm mpi_communicator; + 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 - } + 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; + double + Input::obstacle_function (const double x, + const double y) + { + int ix = (int) (x / hx); + int iy = (int) (y / hy); - if (iy < 0) - iy = 0; + if (ix < 0) + ix = 0; - if (ix >= nx - 1) - ix = nx - 2; + if (iy < 0) + iy = 0; - if (iy >= ny - 1) - iy = ny - 2; + if (ix >= nx - 1) + ix = nx - 2; - double val = 0.0; - { - FullMatrix H(4, 4); - Vector X(4); - Vector b(4); - - double xx = 0.0; - double yy = 0.0; - - 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); - } + if (iy >= ny - 1) + iy = ny - 2; - return val; + double val = 0.0; + { + FullMatrix H(4, 4); + Vector X(4); + Vector b(4); + + double xx = 0.0; + double yy = 0.0; + + 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); } + 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) - { - std::ifstream f(name.c_str()); + void + Input::read_obstacle (const std::string name) + { + std::ifstream f(name.c_str()); - std::string temp; - f >> temp >> nx >> ny; - assert(nx > 0 && ny > 0); + std::string temp; + f >> temp >> nx >> ny; + assert(nx > 0 && ny > 0); - for (int k = 0; k < nx * ny; k++) - { - double val; - f >> val; - obstacle_data.push_back(val); - } + for (int k = 0; k < nx * ny; k++) + { + double val; + f >> val; + obstacle_data.push_back(val); + } - hx = 1.0 / (nx - 1); - hy = 1.0 / (ny - 1); + hx = 1.0 / (nx - 1); + hy = 1.0 / (ny - 1); - pcout << "Resolution of the scanned obstacle picture: " << nx << " x " + pcout << "Resolution of the scanned obstacle picture: " << nx << " x " << ny << std::endl; - } + } // @sect3{The ConstitutiveLaw class template} @@ -268,47 +268,47 @@ namespace Step42 // For gamma = 0 we obtain perfect elastoplastic // behavior. template - class ConstitutiveLaw + class ConstitutiveLaw + { + public: + ConstitutiveLaw ( + double _E, double _nu, double _sigma_0, double _gamma, + MPI_Comm _mpi_communicator, ConditionalOStream _pcout); + + void + plast_linear_hardening ( + SymmetricTensor<4, dim> &stress_strain_tensor, + const SymmetricTensor<2, dim> &strain_tensor, + unsigned int &elast_points, unsigned int &plast_points, + double &yield); + void + linearized_plast_linear_hardening ( + SymmetricTensor<4, dim> &stress_strain_tensor_linearized, + SymmetricTensor<4, dim> &stress_strain_tensor, + const SymmetricTensor<2, dim> &strain_tensor); + inline SymmetricTensor<2, dim> + get_strain ( + const FEValues &fe_values, const unsigned int shape_func, + const unsigned int q_point) const; + void + set_sigma_0 ( + double sigma_hlp) { - public: - ConstitutiveLaw ( - double _E, double _nu, double _sigma_0, double _gamma, - MPI_Comm _mpi_communicator, ConditionalOStream _pcout); - - void - plast_linear_hardening ( - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor, - unsigned int &elast_points, unsigned int &plast_points, - double &yield); - void - linearized_plast_linear_hardening ( - SymmetricTensor<4, dim> &stress_strain_tensor_linearized, - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor); - inline SymmetricTensor<2, dim> - get_strain ( - const FEValues &fe_values, const unsigned int shape_func, - const unsigned int q_point) const; - void - set_sigma_0 ( - double sigma_hlp) - { - sigma_0 = sigma_hlp; - } + sigma_0 = sigma_hlp; + } - private: - SymmetricTensor<4, dim> stress_strain_tensor_mu; - SymmetricTensor<4, dim> stress_strain_tensor_kappa; - double E; - double nu; - double sigma_0; - double gamma; - double mu; - double kappa; - MPI_Comm mpi_communicator; - ConditionalOStream pcout; - }; + private: + SymmetricTensor<4, dim> stress_strain_tensor_mu; + SymmetricTensor<4, dim> stress_strain_tensor_kappa; + double E; + double nu; + double sigma_0; + double gamma; + double mu; + double kappa; + MPI_Comm mpi_communicator; + ConditionalOStream pcout; + }; // The constructor of the ConstitutiveLaw class sets the // required material parameter for our deformable body: @@ -320,45 +320,45 @@ namespace Step42 // of the volumetric and deviator part. For further details // see the documentation above. template - ConstitutiveLaw::ConstitutiveLaw ( - double _E, double _nu, double _sigma_0, double _gamma, - MPI_Comm _mpi_communicator, ConditionalOStream _pcout) - : - E(_E), - nu(_nu), - sigma_0(_sigma_0), - gamma(_gamma), - mpi_communicator(_mpi_communicator), - pcout(_pcout) - { - mu = E / (2 * (1 + nu)); - kappa = E / (3 * (1 - 2 * nu)); - stress_strain_tensor_kappa = kappa - * outer_product(unit_symmetric_tensor(), - unit_symmetric_tensor()); - stress_strain_tensor_mu = 2 * mu - * (identity_tensor() - - outer_product(unit_symmetric_tensor(), - unit_symmetric_tensor()) / 3.0); - } + ConstitutiveLaw::ConstitutiveLaw ( + double _E, double _nu, double _sigma_0, double _gamma, + MPI_Comm _mpi_communicator, ConditionalOStream _pcout) + : + E(_E), + nu(_nu), + sigma_0(_sigma_0), + gamma(_gamma), + mpi_communicator(_mpi_communicator), + pcout(_pcout) + { + mu = E / (2 * (1 + nu)); + kappa = E / (3 * (1 - 2 * nu)); + stress_strain_tensor_kappa = kappa + * outer_product(unit_symmetric_tensor(), + unit_symmetric_tensor()); + stress_strain_tensor_mu = 2 * mu + * (identity_tensor() + - outer_product(unit_symmetric_tensor(), + unit_symmetric_tensor()) / 3.0); + } // @sect3{ConstitutiveLaw::ConstitutiveLaw} // Calculates the strain $\varepsilon(\varphi)=\dfrac{1}{2}\left(\nabla\varphi + \nabla\varphi^T$ // for the shape functions $\varphi$. template - inline 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); - SymmetricTensor<2, dim> tmp; + inline 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); + SymmetricTensor<2, dim> tmp; - tmp = fe_values[displacement].symmetric_gradient(shape_func, q_point); + tmp = fe_values[displacement].symmetric_gradient(shape_func, q_point); - return tmp; - } + return tmp; + } // @sect3{ConstitutiveLaw::plast_linear_hardening} @@ -370,39 +370,39 @@ namespace Step42 // residual in // PlasticityContactProblem::residual_nl_system(TrilinosWrappers::MPI::Vector &u). template - void - ConstitutiveLaw::plast_linear_hardening ( - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor, - unsigned int &elast_points, unsigned int &plast_points, double &yield) - { - if (dim == 3) - { - SymmetricTensor<2, dim> stress_tensor; - stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) - * strain_tensor; + void + ConstitutiveLaw::plast_linear_hardening ( + SymmetricTensor<4, dim> &stress_strain_tensor, + const SymmetricTensor<2, dim> &strain_tensor, + unsigned int &elast_points, unsigned int &plast_points, double &yield) + { + if (dim == 3) + { + SymmetricTensor<2, dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) + * strain_tensor; - SymmetricTensor<2, dim> deviator_stress_tensor = deviator( - stress_tensor); + SymmetricTensor<2, dim> deviator_stress_tensor = deviator( + stress_tensor); - double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); + double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); - yield = 0; - stress_strain_tensor = stress_strain_tensor_mu; - double beta = 1.0; - if (deviator_stress_tensor_norm > sigma_0) - { - beta = sigma_0 / deviator_stress_tensor_norm; - stress_strain_tensor *= (gamma + (1 - gamma) * beta); - yield = 1; - plast_points += 1; - } - else - elast_points += 1; + yield = 0; + stress_strain_tensor = stress_strain_tensor_mu; + double beta = 1.0; + if (deviator_stress_tensor_norm > sigma_0) + { + beta = sigma_0 / deviator_stress_tensor_norm; + stress_strain_tensor *= (gamma + (1 - gamma) * beta); + yield = 1; + plast_points += 1; + } + else + elast_points += 1; - stress_strain_tensor += stress_strain_tensor_kappa; - } - } + stress_strain_tensor += stress_strain_tensor_kappa; + } + } // @sect3{ConstitutiveLaw::linearized_plast_linear_hardening} @@ -417,210 +417,210 @@ namespace Step42 // 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) - { - if (dim == 3) - { - SymmetricTensor<2, dim> stress_tensor; - stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) - * strain_tensor; + 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) + { + if (dim == 3) + { + SymmetricTensor<2, dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) + * strain_tensor; - SymmetricTensor<2, dim> deviator_stress_tensor = deviator( - stress_tensor); + SymmetricTensor<2, dim> deviator_stress_tensor = deviator( + stress_tensor); - double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); + double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); - stress_strain_tensor = stress_strain_tensor_mu; - stress_strain_tensor_linearized = stress_strain_tensor_mu; - double beta = 1.0; - if (deviator_stress_tensor_norm > sigma_0) - { - beta = sigma_0 / deviator_stress_tensor_norm; - stress_strain_tensor *= (gamma + (1 - gamma) * beta); - stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta); - deviator_stress_tensor /= deviator_stress_tensor_norm; - stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu - * outer_product(deviator_stress_tensor, - deviator_stress_tensor); - } + stress_strain_tensor = stress_strain_tensor_mu; + stress_strain_tensor_linearized = stress_strain_tensor_mu; + double beta = 1.0; + if (deviator_stress_tensor_norm > sigma_0) + { + beta = sigma_0 / deviator_stress_tensor_norm; + stress_strain_tensor *= (gamma + (1 - gamma) * beta); + stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta); + deviator_stress_tensor /= deviator_stress_tensor_norm; + stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu + * outer_product(deviator_stress_tensor, + deviator_stress_tensor); + } - stress_strain_tensor += stress_strain_tensor_kappa; - stress_strain_tensor_linearized += stress_strain_tensor_kappa; - } - } + stress_strain_tensor += stress_strain_tensor_kappa; + stress_strain_tensor_linearized += stress_strain_tensor_kappa; + } + } namespace EquationData { // It possible to apply an additional body force // but in here it is set to zero. template - class RightHandSide : public Function + class RightHandSide : public Function + { + public: + RightHandSide () + : + Function(dim) { - public: - RightHandSide () - : - Function(dim) - { - } + } - 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; + double + RightHandSide::value ( + const Point &p, const unsigned int component) const + { + double return_value = 0.0; - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - return_value = 0.0; + if (component == 0) + return_value = 0.0; + if (component == 1) + return_value = 0.0; + if (component == 2) + return_value = 0.0; - return return_value; - } + return return_value; + } template - void - 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); - } + void + 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 + class BoundaryValues : public Function + { + public: + BoundaryValues () + : + Function(dim) { - public: - BoundaryValues () - : - Function(dim) - { - } - ; + } + ; - 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; + double + BoundaryValues::value ( + const Point &p, const unsigned int component) const + { + double return_value = 0; - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - return_value = 0.0; + if (component == 0) + return_value = 0.0; + if (component == 1) + return_value = 0.0; + if (component == 2) + return_value = 0.0; - return return_value; - } + return return_value; + } template - void - 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); - } + void + 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); + } // This function is obviously implemented to // define the obstacle that penetrates our deformable // body. You can choose between two ways to define // your obstacle: to read it from a file or to use // a function (here a ball). -// z_max_domain is the z value of the surface of the work piece +// z_max_domain is the z value of the surface of the work piece template - class Obstacle : public Function + class Obstacle : public Function + { + public: + Obstacle ( + std_cxx1x::shared_ptr > const &_input, + bool _use_read_obstacle, double z_max_domain) + : + Function(dim), + input_obstacle_copy(_input), + use_read_obstacle(_use_read_obstacle), + z_max_domain(z_max_domain) { - public: - Obstacle ( - std_cxx1x::shared_ptr > const &_input, - bool _use_read_obstacle, double z_max_domain) - : - Function(dim), - input_obstacle_copy(_input), - use_read_obstacle(_use_read_obstacle), - z_max_domain(z_max_domain) - { - } + } - 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; - private: - std_cxx1x::shared_ptr > const &input_obstacle_copy; - bool use_read_obstacle; - double z_max_domain; - }; + private: + std_cxx1x::shared_ptr > const &input_obstacle_copy; + bool use_read_obstacle; + double z_max_domain; + }; template - double - Obstacle::value ( - const Point &p, const unsigned int component) const - { - if (component == 0) - return p(0); - if (component == 1) - return p(1); + double + Obstacle::value ( + const Point &p, const unsigned int component) const + { + if (component == 0) + return p(0); + if (component == 1) + return p(1); - //component==2: - if (use_read_obstacle) - { - 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)); - else - return 10000.0; - } - else - { - //sphere: - return -std::sqrt( - 0.36 - (p(0) - 0.5) * (p(0) - 0.5) - - (p(1) - 0.5) * (p(1) - 0.5)) + z_max_domain + 0.59; - } - } + //component==2: + if (use_read_obstacle) + { + 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)); + else + return 10000.0; + } + else + { + //sphere: + return -std::sqrt( + 0.36 - (p(0) - 0.5) * (p(0) - 0.5) + - (p(1) - 0.5) * (p(1) - 0.5)) + z_max_domain + 0.59; + } + } template - void - Obstacle::vector_value ( - const Point &p, Vector &values) const - { - for (unsigned int c = 0; c < this->n_components; ++c) - values(c) = Obstacle::value(p, c); - } + void + Obstacle::vector_value ( + const Point &p, Vector &values) const + { + for (unsigned int c = 0; c < this->n_components; ++c) + values(c) = Obstacle::value(p, c); + } } // @sect3{The PlasticityContactProblem class template} @@ -641,123 +641,123 @@ namespace Step42 // operator for the constitutive law. template - class PlasticityContactProblem + class PlasticityContactProblem + { + public: + PlasticityContactProblem ( + const ParameterHandler &prm); + void + run (); + + static void + declare ( + ParameterHandler &prm); + + private: + void + make_grid (); + void + setup_system (); + void + assemble_nl_system ( + TrilinosWrappers::MPI::Vector &u); + void + residual_nl_system ( + TrilinosWrappers::MPI::Vector &u); + void + assemble_mass_matrix_diagonal ( + TrilinosWrappers::SparseMatrix &mass_matrix); + void + update_solution_and_constraints (); + void + dirichlet_constraints (); + void + solve (); + void + solve_newton (); + void + refine_grid (); + void + move_mesh ( + const TrilinosWrappers::MPI::Vector &_complete_displacement) const; + void + output_results ( + const std::string &title); + void + output_contact_force ( + const unsigned int cycle); + + double to_refine_factor; + double to_coarsen_factor; + unsigned int cycle; + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + FE_Q u; + FESystem fe; + DoFHandler dof_handler; + + // We are using the SolutionTransfer class to interpolate the + // solution on the new refined mesh. It appears in th refine_grid() + // and the run() function. + std_cxx1x::shared_ptr< + parallel::distributed::SolutionTransfer > soltrans; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + unsigned int number_iterations; + + ConstraintMatrix constraints; + ConstraintMatrix constraints_hanging_nodes; + ConstraintMatrix constraints_dirichlet_hanging_nodes; + + TrilinosWrappers::SparseMatrix system_matrix_newton; + + TrilinosWrappers::MPI::Vector solution; + TrilinosWrappers::MPI::Vector system_rhs_newton; + TrilinosWrappers::MPI::Vector system_rhs_lambda; + TrilinosWrappers::MPI::Vector resid_vector; + TrilinosWrappers::MPI::Vector diag_mass_matrix_vector; + Vector cell_constitution; + IndexSet active_set; + + ConditionalOStream pcout; + + TrilinosWrappers::PreconditionAMG::AdditionalData additional_data; + TrilinosWrappers::PreconditionAMG preconditioner_u; + + std_cxx1x::shared_ptr > input_obstacle; + std_cxx1x::shared_ptr > plast_lin_hard; + + double sigma_0; // Yield stress + double gamma; // Parameter for the linear isotropic hardening + double e_modul; // E-Modul + double nu; // Poisson ratio + + TimerOutput computing_timer; + + unsigned int degree; + unsigned int n_initial_refinements; + struct RefinementStrategy { - public: - PlasticityContactProblem ( - const ParameterHandler &prm); - void - run (); - - static void - declare ( - ParameterHandler &prm); - - private: - void - make_grid (); - void - setup_system (); - void - assemble_nl_system ( - TrilinosWrappers::MPI::Vector &u); - void - residual_nl_system ( - TrilinosWrappers::MPI::Vector &u); - void - assemble_mass_matrix_diagonal ( - TrilinosWrappers::SparseMatrix &mass_matrix); - void - update_solution_and_constraints (); - void - dirichlet_constraints (); - void - solve (); - void - solve_newton (); - void - refine_grid (); - void - move_mesh ( - const TrilinosWrappers::MPI::Vector &_complete_displacement) const; - void - output_results ( - const std::string &title); - void - output_contact_force ( - const unsigned int cycle); - - double to_refine_factor; - double to_coarsen_factor; - unsigned int cycle; - - MPI_Comm mpi_communicator; - - parallel::distributed::Triangulation triangulation; - - FE_Q u; - FESystem fe; - DoFHandler dof_handler; - - // We are using the SolutionTransfer class to interpolate the - // solution on the new refined mesh. It appears in th refine_grid() - // and the run() function. - std_cxx1x::shared_ptr< - parallel::distributed::SolutionTransfer > soltrans; - - IndexSet locally_owned_dofs; - IndexSet locally_relevant_dofs; - - unsigned int number_iterations; - - ConstraintMatrix constraints; - ConstraintMatrix constraints_hanging_nodes; - ConstraintMatrix constraints_dirichlet_hanging_nodes; - - TrilinosWrappers::SparseMatrix system_matrix_newton; - - TrilinosWrappers::MPI::Vector solution; - TrilinosWrappers::MPI::Vector system_rhs_newton; - TrilinosWrappers::MPI::Vector system_rhs_lambda; - TrilinosWrappers::MPI::Vector resid_vector; - TrilinosWrappers::MPI::Vector diag_mass_matrix_vector; - Vector cell_constitution; - IndexSet active_set; - - ConditionalOStream pcout; - - TrilinosWrappers::PreconditionAMG::AdditionalData additional_data; - TrilinosWrappers::PreconditionAMG preconditioner_u; - - std_cxx1x::shared_ptr > input_obstacle; - std_cxx1x::shared_ptr > plast_lin_hard; - - double sigma_0; // Yield stress - double gamma; // Parameter for the linear isotropic hardening - double e_modul; // E-Modul - double nu; // Poisson ratio - - TimerOutput computing_timer; - - unsigned int degree; - unsigned int n_initial_refinements; - struct RefinementStrategy - { - enum value - { - refine_global, - refine_percentage, - refine_fix_dofs - }; - }; - typename RefinementStrategy::value refinement_strategy; - unsigned int n_cycles; - std::string obstacle_filename; - std::string output_dir; - bool transfer_solution; - std::string base_mesh; + enum value + { + refine_global, + refine_percentage, + refine_fix_dofs + }; }; + typename RefinementStrategy::value refinement_strategy; + unsigned int n_cycles; + std::string obstacle_filename; + std::string output_dir; + bool transfer_solution; + std::string base_mesh; + }; // @sect3{Implementation of the PlasticityContactProblem class} @@ -766,90 +766,90 @@ namespace Step42 // above. As before, we will write everything template - PlasticityContactProblem::PlasticityContactProblem ( - const ParameterHandler &prm) - : - mpi_communicator(MPI_COMM_WORLD), - triangulation(mpi_communicator), - u(QGaussLobatto<1>(prm.get_integer("polynomial degree") + 1)), - fe(u, dim), - dof_handler(triangulation), - pcout(std::cout, - (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), - sigma_0(400.0), - gamma(0.01), - e_modul(2.0e+5), - nu(0.3), - computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never, - TimerOutput::wall_times) - { - // double _E, double _nu, double _sigma_0, double _gamma - plast_lin_hard.reset( - new ConstitutiveLaw(e_modul, nu, sigma_0, gamma, - mpi_communicator, pcout)); - - degree = prm.get_integer("polynomial degree"); - n_initial_refinements = prm.get_integer("number of initial refinements"); - std::string strat = prm.get("refinement strategy"); - if (strat == "global") - refinement_strategy = RefinementStrategy::refine_global; - else if (strat == "percentage") - refinement_strategy = RefinementStrategy::refine_percentage; - else - throw ExcNotImplemented(); - - n_cycles = prm.get_integer("number of cycles"); - obstacle_filename = prm.get("obstacle filename"); - output_dir = prm.get("output directory"); - if (output_dir != "" && *(output_dir.rbegin()) != '/') - output_dir += "/"; - mkdir(output_dir.c_str(), 0777); - - transfer_solution = prm.get_bool("transfer solution"); - base_mesh = prm.get("base mesh"); - - pcout << " Using output directory '" << output_dir << "'" << std::endl; - pcout << " FE degree " << degree << std::endl; - pcout << " Obstacle '" << obstacle_filename << "'" << std::endl; - pcout << " transfer solution " + PlasticityContactProblem::PlasticityContactProblem ( + const ParameterHandler &prm) + : + mpi_communicator(MPI_COMM_WORLD), + triangulation(mpi_communicator), + u(QGaussLobatto<1>(prm.get_integer("polynomial degree") + 1)), + fe(u, dim), + dof_handler(triangulation), + pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), + sigma_0(400.0), + gamma(0.01), + e_modul(2.0e+5), + nu(0.3), + computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never, + TimerOutput::wall_times) + { + // double _E, double _nu, double _sigma_0, double _gamma + plast_lin_hard.reset( + new ConstitutiveLaw(e_modul, nu, sigma_0, gamma, + mpi_communicator, pcout)); + + degree = prm.get_integer("polynomial degree"); + n_initial_refinements = prm.get_integer("number of initial refinements"); + std::string strat = prm.get("refinement strategy"); + if (strat == "global") + refinement_strategy = RefinementStrategy::refine_global; + else if (strat == "percentage") + refinement_strategy = RefinementStrategy::refine_percentage; + else + throw ExcNotImplemented(); + + n_cycles = prm.get_integer("number of cycles"); + obstacle_filename = prm.get("obstacle filename"); + output_dir = prm.get("output directory"); + if (output_dir != "" && *(output_dir.rbegin()) != '/') + output_dir += "/"; + mkdir(output_dir.c_str(), 0777); + + transfer_solution = prm.get_bool("transfer solution"); + base_mesh = prm.get("base mesh"); + + pcout << " Using output directory '" << output_dir << "'" << std::endl; + pcout << " FE degree " << degree << std::endl; + pcout << " Obstacle '" << obstacle_filename << "'" << std::endl; + pcout << " transfer solution " << (transfer_solution ? "true" : "false") << std::endl; - } + } // @sect4{PlasticityContactProblem::declare} template - void - PlasticityContactProblem::declare ( - ParameterHandler &prm) - { - prm.declare_entry("polynomial degree", "1", Patterns::Integer(), - "polynomial degree of the FE_Q finite element space, typically 1 or 2"); - prm.declare_entry("number of initial refinements", "2", - Patterns::Integer(), - "number of initial global refinements before the first computation"); - prm.declare_entry("refinement strategy", "percentage", - Patterns::Selection("global|percentage|fix dofs"), - "refinement strategy for each cycle:\n" - " global: one global refinement\n" - "percentage: fixed percentage gets refined using kelly\n" - " fix dofs: tries to achieve 2^initial_refinement*300 dofs after cycle 1 (only use 2 cycles!). Changes the coarse mesh!"); - prm.declare_entry("number of cycles", "5", Patterns::Integer(), - "number of adaptive cycles to run"); - prm.declare_entry("obstacle filename", "", Patterns::Anything(), - "obstacle file to read, use 'obstacle_file.pbm' or leave empty to use a sphere"); - prm.declare_entry("output directory", "", Patterns::Anything(), - "directory to put output files (graphical output and benchmark statistics), leave empty to put into current directory"); - prm.declare_entry("transfer solution", "false", Patterns::Bool(), - "decide if the solution should be used as a starting guess for the finer mesh, use 0 otherwise."); - prm.declare_entry("base mesh", "box", - Patterns::Selection("box|half sphere"), - "select the shape of the work piece: 'box' or 'half sphere'"); + void + PlasticityContactProblem::declare ( + ParameterHandler &prm) + { + prm.declare_entry("polynomial degree", "1", Patterns::Integer(), + "polynomial degree of the FE_Q finite element space, typically 1 or 2"); + prm.declare_entry("number of initial refinements", "2", + Patterns::Integer(), + "number of initial global refinements before the first computation"); + prm.declare_entry("refinement strategy", "percentage", + Patterns::Selection("global|percentage|fix dofs"), + "refinement strategy for each cycle:\n" + " global: one global refinement\n" + "percentage: fixed percentage gets refined using kelly\n" + " fix dofs: tries to achieve 2^initial_refinement*300 dofs after cycle 1 (only use 2 cycles!). Changes the coarse mesh!"); + prm.declare_entry("number of cycles", "5", Patterns::Integer(), + "number of adaptive cycles to run"); + prm.declare_entry("obstacle filename", "", Patterns::Anything(), + "obstacle file to read, use 'obstacle_file.pbm' or leave empty to use a sphere"); + prm.declare_entry("output directory", "", Patterns::Anything(), + "directory to put output files (graphical output and benchmark statistics), leave empty to put into current directory"); + prm.declare_entry("transfer solution", "false", Patterns::Bool(), + "decide if the solution should be used as a starting guess for the finer mesh, use 0 otherwise."); + prm.declare_entry("base mesh", "box", + Patterns::Selection("box|half sphere"), + "select the shape of the work piece: 'box' or 'half sphere'"); - } + } Point<3> rotate_half_sphere ( - const Point<3> &in) + const Point<3> &in) { return Point<3>(in(2), in(1), -in(0)); } @@ -857,600 +857,600 @@ namespace Step42 // @sect4{PlasticityContactProblem::make_grid} template - void - PlasticityContactProblem::make_grid () - { - - if (base_mesh == "half sphere") - { - Point center(0, 0, 0); - double radius = 0.8; - GridGenerator::half_hyper_ball(triangulation, center, radius); - GridTools::transform(&rotate_half_sphere, triangulation); - Point shift(0.5, 0.5, 0.5); - GridTools::shift(shift, triangulation); - static HyperBallBoundary boundary_description( - Point(0.5, 0.5, 0.5), radius); - triangulation.set_boundary(0, boundary_description); - - triangulation.refine_global(n_initial_refinements); - - to_refine_factor = 0.3; - to_coarsen_factor = 0.03; - return; - } - - Point p1(0, 0, 0); - Point p2(1.0, 1.0, 1.0); + void + PlasticityContactProblem::make_grid () + { - GridGenerator::hyper_rectangle(triangulation, p1, p2); - to_refine_factor = 0.3; - to_coarsen_factor = 0.03; + if (base_mesh == "half sphere") + { + Point center(0, 0, 0); + double radius = 0.8; + GridGenerator::half_hyper_ball(triangulation, center, radius); + GridTools::transform(&rotate_half_sphere, triangulation); + Point shift(0.5, 0.5, 0.5); + GridTools::shift(shift, triangulation); + static HyperBallBoundary boundary_description( + Point(0.5, 0.5, 0.5), radius); + triangulation.set_boundary(0, boundary_description); + + triangulation.refine_global(n_initial_refinements); + + to_refine_factor = 0.3; + to_coarsen_factor = 0.03; + return; + } - Triangulation<3>::active_cell_iterator cell = - triangulation.begin_active(), endc = triangulation.end(); + Point p1(0, 0, 0); + Point p2(1.0, 1.0, 1.0); - /* boundary_indicators: - _______ - / 1 /| - /______ / | - | | 8| - | 8 | / - |_______|/ - 6 + GridGenerator::hyper_rectangle(triangulation, p1, p2); + to_refine_factor = 0.3; + to_coarsen_factor = 0.03; - The boundary indicators of the sides of the cube are 8. - The boundary indicator of the bottom is indicated with 6 - and the top with 1. - */ + Triangulation<3>::active_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 (cell->face(face)->center()[2] == p2(2)) - cell->face(face)->set_boundary_indicator(1); - if (cell->face(face)->center()[0] == p1(0) - || cell->face(face)->center()[0] == p2(0) - || cell->face(face)->center()[1] == p1(1) - || cell->face(face)->center()[1] == p2(1)) - cell->face(face)->set_boundary_indicator(8); - if (cell->face(face)->center()[2] == p1(2)) - cell->face(face)->set_boundary_indicator(6); - } + /* boundary_indicators: + _______ + / 1 /| + /______ / | + | | 8| + | 8 | / + |_______|/ + 6 + + The boundary indicators of the sides of the cube are 8. + The boundary indicator of the bottom is indicated with 6 + and the top with 1. + */ + + for (; cell != endc; ++cell) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + if (cell->face(face)->center()[2] == p2(2)) + cell->face(face)->set_boundary_indicator(1); + if (cell->face(face)->center()[0] == p1(0) + || cell->face(face)->center()[0] == p2(0) + || cell->face(face)->center()[1] == p1(1) + || cell->face(face)->center()[1] == p2(1)) + cell->face(face)->set_boundary_indicator(8); + if (cell->face(face)->center()[2] == p1(2)) + cell->face(face)->set_boundary_indicator(6); + } - triangulation.refine_global(n_initial_refinements); - } + triangulation.refine_global(n_initial_refinements); + } template - void - PlasticityContactProblem::setup_system () + void + PlasticityContactProblem::setup_system () + { + // setup dofs { - // setup dofs - { - TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs"); - dof_handler.distribute_dofs(fe); + TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs"); + dof_handler.distribute_dofs(fe); - locally_owned_dofs = dof_handler.locally_owned_dofs(); - locally_relevant_dofs.clear(); - DoFTools::extract_locally_relevant_dofs(dof_handler, - locally_relevant_dofs); - } + locally_owned_dofs = dof_handler.locally_owned_dofs(); + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs(dof_handler, + locally_relevant_dofs); + } - // setup hanging nodes and Dirichlet constraints - { - TimerOutput::Scope t(computing_timer, "Setup: constraints"); - constraints_hanging_nodes.reinit(locally_relevant_dofs); - DoFTools::make_hanging_node_constraints(dof_handler, - constraints_hanging_nodes); - constraints_hanging_nodes.close(); - - pcout << " Number of active cells: " - << triangulation.n_global_active_cells() << std::endl - << " Number of degrees of freedom: " << dof_handler.n_dofs() - << std::endl; - - dirichlet_constraints(); - } + // setup hanging nodes and Dirichlet constraints + { + TimerOutput::Scope t(computing_timer, "Setup: constraints"); + constraints_hanging_nodes.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, + constraints_hanging_nodes); + constraints_hanging_nodes.close(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + dirichlet_constraints(); + } - // Initialization for matrices and vectors - { - TimerOutput::Scope t(computing_timer, "Setup: vectors"); - solution.reinit(locally_relevant_dofs, mpi_communicator); - system_rhs_newton.reinit(locally_owned_dofs, mpi_communicator); - system_rhs_lambda.reinit(system_rhs_newton); - resid_vector.reinit(system_rhs_newton); - diag_mass_matrix_vector.reinit(system_rhs_newton); - cell_constitution.reinit(triangulation.n_active_cells()); - active_set.clear(); - active_set.set_size(locally_relevant_dofs.size()); - } + // Initialization for matrices and vectors + { + TimerOutput::Scope t(computing_timer, "Setup: vectors"); + solution.reinit(locally_relevant_dofs, mpi_communicator); + system_rhs_newton.reinit(locally_owned_dofs, mpi_communicator); + system_rhs_lambda.reinit(system_rhs_newton); + resid_vector.reinit(system_rhs_newton); + diag_mass_matrix_vector.reinit(system_rhs_newton); + cell_constitution.reinit(triangulation.n_active_cells()); + active_set.clear(); + active_set.set_size(locally_relevant_dofs.size()); + } - // setup sparsity pattern - { - TimerOutput::Scope t(computing_timer, "Setup: matrix"); - TrilinosWrappers::SparsityPattern sp(locally_owned_dofs, - mpi_communicator); + // setup sparsity pattern + { + TimerOutput::Scope t(computing_timer, "Setup: matrix"); + TrilinosWrappers::SparsityPattern sp(locally_owned_dofs, + mpi_communicator); - DoFTools::make_sparsity_pattern(dof_handler, sp, - constraints_dirichlet_hanging_nodes, false, - Utilities::MPI::this_mpi_process(mpi_communicator)); + DoFTools::make_sparsity_pattern(dof_handler, sp, + constraints_dirichlet_hanging_nodes, false, + Utilities::MPI::this_mpi_process(mpi_communicator)); - sp.compress(); + sp.compress(); - system_matrix_newton.reinit(sp); + system_matrix_newton.reinit(sp); - // we are going to reuse the system - // matrix for assembling the diagonal - // of the mass matrix so that we do not - // need to allocate two sparse matrices - // at the same time: - TrilinosWrappers::SparseMatrix & mass_matrix = system_matrix_newton; - assemble_mass_matrix_diagonal(mass_matrix); - const unsigned int start = (system_rhs_newton.local_range().first), - end = (system_rhs_newton.local_range().second); - for (unsigned int j = start; j < end; j++) - diag_mass_matrix_vector(j) = mass_matrix.diag_element(j); + // we are going to reuse the system + // matrix for assembling the diagonal + // of the mass matrix so that we do not + // need to allocate two sparse matrices + // at the same time: + TrilinosWrappers::SparseMatrix &mass_matrix = system_matrix_newton; + assemble_mass_matrix_diagonal(mass_matrix); + const unsigned int start = (system_rhs_newton.local_range().first), + end = (system_rhs_newton.local_range().second); + for (unsigned int j = start; j < end; j++) + diag_mass_matrix_vector(j) = mass_matrix.diag_element(j); - number_iterations = 0; + number_iterations = 0; - diag_mass_matrix_vector.compress(VectorOperation::insert); + diag_mass_matrix_vector.compress(VectorOperation::insert); - // remove the mass matrix entries from the matrix: - mass_matrix = 0; - } + // remove the mass matrix entries from the matrix: + mass_matrix = 0; } + } template - void - PlasticityContactProblem::assemble_nl_system ( - TrilinosWrappers::MPI::Vector &u) - { - TimerOutput::Scope t(computing_timer, "Assembling"); + void + PlasticityContactProblem::assemble_nl_system ( + TrilinosWrappers::MPI::Vector &u) + { + TimerOutput::Scope t(computing_timer, "Assembling"); - QGauss quadrature_formula(fe.degree + 1); - QGauss face_quadrature_formula(fe.degree + 1); + QGauss quadrature_formula(fe.degree + 1); + QGauss face_quadrature_formula(fe.degree + 1); - FEValues fe_values(fe, quadrature_formula, - UpdateFlags( - update_values | update_gradients | update_q_points - | update_JxW_values)); + FEValues fe_values(fe, quadrature_formula, + UpdateFlags( + update_values | update_gradients | update_q_points + | update_JxW_values)); - FEFaceValues fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points | update_JxW_values); + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); - const EquationData::RightHandSide right_hand_side; - std::vector > right_hand_side_values(n_q_points, - Vector(dim)); - std::vector > right_hand_side_values_face(n_face_q_points, - Vector(dim)); + const EquationData::RightHandSide right_hand_side; + std::vector > right_hand_side_values(n_q_points, + Vector(dim)); + std::vector > right_hand_side_values_face(n_face_q_points, + Vector(dim)); - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - Vector cell_rhs(dofs_per_cell); + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); - std::vector local_dof_indices(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); - const FEValuesExtractors::Vector displacement(0); + const FEValuesExtractors::Vector displacement(0); - const double kappa = 1.0; - for (; cell != endc; ++cell) - if (cell->is_locally_owned()) - { - fe_values.reinit(cell); - cell_matrix = 0; - cell_rhs = 0; + const double kappa = 1.0; + for (; cell != endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - right_hand_side_values); + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + right_hand_side_values); - std::vector > strain_tensor(n_q_points); - fe_values[displacement].get_function_symmetric_gradients(u, - strain_tensor); + std::vector > strain_tensor(n_q_points); + fe_values[displacement].get_function_symmetric_gradients(u, + strain_tensor); - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - { - SymmetricTensor<4, dim> stress_strain_tensor_linearized; - SymmetricTensor<4, dim> stress_strain_tensor; - SymmetricTensor<2, dim> stress_tensor; + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + SymmetricTensor<4, dim> stress_strain_tensor_linearized; + SymmetricTensor<4, dim> stress_strain_tensor; + 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->linearized_plast_linear_hardening( + stress_strain_tensor_linearized, stress_strain_tensor, + strain_tensor[q_point]); - 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); - - 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.JxW(q_point)); - } - - // the linearized part a(v^i;v^i,v) of the rhs - cell_rhs(i) += (stress_tensor * strain_tensor[q_point] - * fe_values.JxW(q_point)); - - // 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.JxW(q_point)); - - // the residual part F(v) of the rhs - Tensor<1, dim> rhs_values; - rhs_values = 0; - cell_rhs(i) += (fe_values[displacement].value(i, q_point) - * rhs_values * fe_values.JxW(q_point)); - } - } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + stress_tensor = stress_strain_tensor_linearized + * plast_lin_hard->get_strain(fe_values, 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.JxW(q_point)); + } + + // the linearized part a(v^i;v^i,v) of the rhs + cell_rhs(i) += (stress_tensor * strain_tensor[q_point] + * fe_values.JxW(q_point)); + + // 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.JxW(q_point)); + + // the residual part F(v) of the rhs + Tensor<1, dim> rhs_values; + rhs_values = 0; + cell_rhs(i) += (fe_values[displacement].value(i, q_point) + * rhs_values * fe_values.JxW(q_point)); + } + } + + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && cell->face(face)->boundary_indicator() == 1) + { + fe_values_face.reinit(cell, face); + + right_hand_side.vector_value_list( + fe_values_face.get_quadrature_points(), + right_hand_side_values_face); + + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) + { + Tensor<1, dim> rhs_values; + rhs_values[2] = right_hand_side_values[q_point][2]; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + cell_rhs(i) += (fe_values_face[displacement].value(i, + q_point) * rhs_values + * fe_values_face.JxW(q_point)); + } + } + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, cell_rhs, + local_dof_indices, system_matrix_newton, system_rhs_newton, + true); + + }; + + system_matrix_newton.compress(VectorOperation::add); + system_rhs_newton.compress(VectorOperation::add); + } + + template + void + PlasticityContactProblem::residual_nl_system ( + TrilinosWrappers::MPI::Vector &u) + { + QGauss quadrature_formula(fe.degree + 1); + QGauss face_quadrature_formula(fe.degree + 1); + + FEValues fe_values(fe, quadrature_formula, + UpdateFlags( + update_values | update_gradients | update_q_points + | update_JxW_values)); - for (unsigned int face = 0; - face < GeometryInfo::faces_per_cell; ++face) - { - if (cell->face(face)->at_boundary() - && cell->face(face)->boundary_indicator() == 1) - { - fe_values_face.reinit(cell, face); - - right_hand_side.vector_value_list( - fe_values_face.get_quadrature_points(), - right_hand_side_values_face); - - for (unsigned int q_point = 0; q_point < n_face_q_points; - ++q_point) - { - Tensor<1, dim> rhs_values; - rhs_values[2] = right_hand_side_values[q_point][2]; - 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)); - } - } - } + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); - cell->get_dof_indices(local_dof_indices); - constraints.distribute_local_to_global(cell_matrix, cell_rhs, - local_dof_indices, system_matrix_newton, system_rhs_newton, - true); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); - }; + const EquationData::RightHandSide right_hand_side; + std::vector > right_hand_side_values(n_q_points, + Vector(dim)); + std::vector > right_hand_side_values_face(n_face_q_points, + Vector(dim)); - system_matrix_newton.compress(VectorOperation::add); - system_rhs_newton.compress(VectorOperation::add); - } + Vector cell_rhs(dofs_per_cell); - template - void - PlasticityContactProblem::residual_nl_system ( - TrilinosWrappers::MPI::Vector &u) - { - QGauss quadrature_formula(fe.degree + 1); - QGauss face_quadrature_formula(fe.degree + 1); + std::vector local_dof_indices(dofs_per_cell); - FEValues fe_values(fe, quadrature_formula, - UpdateFlags( - update_values | update_gradients | update_q_points - | update_JxW_values)); + const FEValuesExtractors::Vector displacement(0); - FEFaceValues fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points | update_JxW_values); + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); + unsigned int elast_points = 0; + unsigned int plast_points = 0; + double yield = 0; + unsigned int cell_number = 0; + cell_constitution = 0; - const EquationData::RightHandSide right_hand_side; - std::vector > right_hand_side_values(n_q_points, - Vector(dim)); - std::vector > right_hand_side_values_face(n_face_q_points, - Vector(dim)); + for (; cell != endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_rhs = 0; - Vector cell_rhs(dofs_per_cell); + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + right_hand_side_values); - std::vector local_dof_indices(dofs_per_cell); + std::vector > strain_tensor(n_q_points); + fe_values[displacement].get_function_symmetric_gradients(u, + strain_tensor); - const FEValuesExtractors::Vector displacement(0); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + SymmetricTensor<4, dim> stress_strain_tensor; + SymmetricTensor<2, dim> stress_tensor; - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); + plast_lin_hard->plast_linear_hardening(stress_strain_tensor, + strain_tensor[q_point], elast_points, plast_points, yield); - unsigned int elast_points = 0; - unsigned int plast_points = 0; - double yield = 0; - unsigned int cell_number = 0; - cell_constitution = 0; + cell_constitution(cell_number) += yield; + 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.JxW(q_point)); + + Tensor<1, dim> rhs_values; + rhs_values = 0; + cell_rhs(i) += ((fe_values[displacement].value(i, q_point) + * rhs_values) * fe_values.JxW(q_point)); + }; + }; - for (; cell != endc; ++cell) - if (cell->is_locally_owned()) - { - fe_values.reinit(cell); - cell_rhs = 0; + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() + && cell->face(face)->boundary_indicator() == 1) + { + fe_values_face.reinit(cell, face); - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - right_hand_side_values); + right_hand_side.vector_value_list( + fe_values_face.get_quadrature_points(), + right_hand_side_values_face); - std::vector > strain_tensor(n_q_points); - fe_values[displacement].get_function_symmetric_gradients(u, - strain_tensor); + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) + { + Tensor<1, dim> rhs_values; + rhs_values[2] = right_hand_side_values[q_point][2]; + 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)); + } + } + } - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - { - SymmetricTensor<4, dim> stress_strain_tensor; - SymmetricTensor<2, dim> stress_tensor; + cell->get_dof_indices(local_dof_indices); + constraints_dirichlet_hanging_nodes.distribute_local_to_global( + cell_rhs, local_dof_indices, system_rhs_newton); - plast_lin_hard->plast_linear_hardening(stress_strain_tensor, - strain_tensor[q_point], elast_points, plast_points, yield); + for (unsigned int i = 0; i < dofs_per_cell; i++) + system_rhs_lambda(local_dof_indices[i]) += cell_rhs(i); - cell_constitution(cell_number) += yield; - 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.JxW(q_point)); - - Tensor<1, dim> rhs_values; - rhs_values = 0; - cell_rhs(i) += ((fe_values[displacement].value(i, q_point) - * rhs_values) * fe_values.JxW(q_point)); - }; - }; - - for (unsigned int face = 0; - face < GeometryInfo::faces_per_cell; ++face) - { - if (cell->face(face)->at_boundary() - && cell->face(face)->boundary_indicator() == 1) - { - fe_values_face.reinit(cell, face); - - right_hand_side.vector_value_list( - fe_values_face.get_quadrature_points(), - right_hand_side_values_face); - - for (unsigned int q_point = 0; q_point < n_face_q_points; - ++q_point) - { - Tensor<1, dim> rhs_values; - rhs_values[2] = right_hand_side_values[q_point][2]; - for (unsigned int i = 0; i < dofs_per_cell; ++i) - cell_rhs(i) += (fe_values_face[displacement].value(i, - q_point) * rhs_values - * fe_values_face.JxW(q_point)); - } - } - } + cell_number += 1; + } + else + { + cell_constitution(cell_number) = 0; + cell_number += 1; + }; - cell->get_dof_indices(local_dof_indices); - constraints_dirichlet_hanging_nodes.distribute_local_to_global( - cell_rhs, local_dof_indices, system_rhs_newton); + cell_constitution /= n_q_points; + cell_constitution.compress(VectorOperation::add); + system_rhs_newton.compress(VectorOperation::add); + system_rhs_lambda.compress(VectorOperation::add); - for (unsigned int i = 0; i < dofs_per_cell; i++) - system_rhs_lambda(local_dof_indices[i]) += cell_rhs(i); +// constraints_hanging_nodes.condense(system_rhs_lambda); - cell_number += 1; - } - else - { - cell_constitution(cell_number) = 0; - cell_number += 1; - }; - - cell_constitution /= n_q_points; - cell_constitution.compress(VectorOperation::add); - system_rhs_newton.compress(VectorOperation::add); - system_rhs_lambda.compress(VectorOperation::add); - -// constraints_hanging_nodes.condense(system_rhs_lambda); - - unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, - mpi_communicator); - unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, - mpi_communicator); - pcout << " Number of elastic quadrature points: " << sum_elast_points + unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, + mpi_communicator); + unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, + mpi_communicator); + pcout << " Number of elastic quadrature points: " << sum_elast_points << " and plastic quadrature points: " << sum_plast_points << std::endl; - } + } template - void - PlasticityContactProblem::assemble_mass_matrix_diagonal ( - TrilinosWrappers::SparseMatrix &mass_matrix) - { - QGaussLobatto face_quadrature_formula(fe.degree + 1); + void + PlasticityContactProblem::assemble_mass_matrix_diagonal ( + TrilinosWrappers::SparseMatrix &mass_matrix) + { + QGaussLobatto face_quadrature_formula(fe.degree + 1); - FEFaceValues fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points | update_JxW_values); + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_face_q_points = face_quadrature_formula.size(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_face_q_points = face_quadrature_formula.size(); - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - Tensor<1, dim, double> ones(dim); - for (unsigned i = 0; i < dim; i++) - ones[i] = 1.0; + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Tensor<1, dim, double> ones(dim); + for (unsigned i = 0; i < dim; i++) + ones[i] = 1.0; - std::vector local_dof_indices(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); - const FEValuesExtractors::Vector displacement(0); + const FEValuesExtractors::Vector displacement(0); - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); + 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 = 0; face < GeometryInfo::faces_per_cell; - ++face) - if (cell->face(face)->at_boundary() - && cell->face(face)->boundary_indicator() == 1) - { - fe_values_face.reinit(cell, face); - cell_matrix = 0; + for (; cell != endc; ++cell) + if (cell->is_locally_owned()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + if (cell->face(face)->at_boundary() + && cell->face(face)->boundary_indicator() == 1) + { + fe_values_face.reinit(cell, face); + cell_matrix = 0; - for (unsigned int q_point = 0; q_point < n_face_q_points; - ++q_point) - for (unsigned int i = 0; i < dofs_per_cell; ++i) - cell_matrix(i, i) += (fe_values_face[displacement].value(i, - q_point) * ones * fe_values_face.JxW(q_point)); + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + cell_matrix(i, i) += (fe_values_face[displacement].value(i, + q_point) * ones * fe_values_face.JxW(q_point)); - cell->get_dof_indices(local_dof_indices); + cell->get_dof_indices(local_dof_indices); -// constraints_dirichlet_hanging_nodes.distribute_local_to_global( -// cell_matrix, local_dof_indices, mass_matrix); +// constraints_dirichlet_hanging_nodes.distribute_local_to_global( +// cell_matrix, local_dof_indices, mass_matrix); - for (unsigned int i = 0; i < dofs_per_cell; i++) - mass_matrix.add(local_dof_indices[i], local_dof_indices[i], - cell_matrix(i, i)); - } - mass_matrix.compress(VectorOperation::add); - } + for (unsigned int i = 0; i < dofs_per_cell; i++) + mass_matrix.add(local_dof_indices[i], local_dof_indices[i], + cell_matrix(i, i)); + } + mass_matrix.compress(VectorOperation::add); + } // @sect4{PlasticityContactProblem::update_solution_and_constraints} // Projection and updating of the active set // for the dofs which penetrates the obstacle. template - void - PlasticityContactProblem::update_solution_and_constraints () - { - const EquationData::Obstacle obstacle(input_obstacle, - (obstacle_filename != ""), (base_mesh == "box" ? 1.0 : 0.5)); - std::vector vertex_touched(dof_handler.n_dofs(), false); - - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); - - TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton); - distributed_solution = solution; - TrilinosWrappers::MPI::Vector lambda(solution); - lambda = resid_vector; - TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(solution); - diag_mass_matrix_vector_relevant = diag_mass_matrix_vector; - - constraints.reinit(locally_relevant_dofs); - active_set.clear(); - IndexSet active_set_locally_owned; - active_set_locally_owned.set_size(locally_owned_dofs.size()); - const double c = 100.0 * e_modul; - - Quadrature face_quadrature(fe.get_unit_face_support_points()); - FEFaceValues fe_values_face(fe, face_quadrature, - update_quadrature_points); - - const unsigned int dofs_per_face = fe.dofs_per_face; - const unsigned int n_face_q_points = face_quadrature.size(); + void + PlasticityContactProblem::update_solution_and_constraints () + { + const EquationData::Obstacle obstacle(input_obstacle, + (obstacle_filename != ""), (base_mesh == "box" ? 1.0 : 0.5)); + std::vector vertex_touched(dof_handler.n_dofs(), false); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); + + TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton); + distributed_solution = solution; + TrilinosWrappers::MPI::Vector lambda(solution); + lambda = resid_vector; + TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(solution); + diag_mass_matrix_vector_relevant = diag_mass_matrix_vector; + + constraints.reinit(locally_relevant_dofs); + active_set.clear(); + IndexSet active_set_locally_owned; + active_set_locally_owned.set_size(locally_owned_dofs.size()); + const double c = 100.0 * e_modul; + + Quadrature face_quadrature(fe.get_unit_face_support_points()); + FEFaceValues fe_values_face(fe, face_quadrature, + update_quadrature_points); + + const unsigned int dofs_per_face = fe.dofs_per_face; + const unsigned int n_face_q_points = face_quadrature.size(); + + // pcout<< "dofs_per_face = " << dofs_per_face + // << "n_face_q_points = " << n_face_q_points + // <is_artificial()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + if (cell->face(face)->at_boundary() + && cell->face(face)->boundary_indicator() == 1) + { + fe_values_face.reinit(cell, face); + std::vector dof_indices(dofs_per_face); + cell->face(face)->get_dof_indices(dof_indices); - // pcout<< "dofs_per_face = " << dofs_per_face - // << "n_face_q_points = " << n_face_q_points - // <is_artificial()) - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; - ++face) - if (cell->face(face)->at_boundary() - && cell->face(face)->boundary_indicator() == 1) - { - fe_values_face.reinit(cell, face); - std::vector dof_indices(dofs_per_face); - cell->face(face)->get_dof_indices(dof_indices); + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) + { + unsigned int component = fe.face_system_to_component_index( + q_point).first; - for (unsigned int q_point = 0; q_point < n_face_q_points; - ++q_point) - { - unsigned int component = fe.face_system_to_component_index( - q_point).first; - - if (component == 2) - { - unsigned int index_z = dof_indices[q_point]; - - if (vertex_touched[index_z] == false) - vertex_touched[index_z] = true; - else - continue; - - // the local row where - Point point( - fe_values_face.quadrature_point(q_point)); - - double obstacle_value = obstacle.value(point, 2); - double solution_index_z = solution(index_z); - double gap = obstacle_value - point(2); - - if (lambda(index_z) - / diag_mass_matrix_vector_relevant(index_z) - + c * (solution_index_z - gap) > 0 - && !(constraints_hanging_nodes.is_constrained( - index_z))) - { - constraints.add_line(index_z); - constraints.set_inhomogeneity(index_z, gap); - distributed_solution(index_z) = gap; - - if (locally_owned_dofs.is_element(index_z)) - { - active_set_locally_owned.add_index(index_z); - if (locally_relevant_dofs.is_element(index_z)) - active_set.add_index(index_z); - } - - } - else if (lambda(index_z) - / diag_mass_matrix_vector_relevant(index_z) - + c * (solution_index_z - gap) > 0 - && constraints_hanging_nodes.is_constrained( - index_z)) - { - if (locally_owned_dofs.is_element(index_z)) - { - counter_hanging_nodes += 1; - -// std::cout << "index_z = " << index_z -// << ", lambda = " << lambda (index_z) -// << ", solution_index_z - gap = " << solution_index_z - gap -// << ", diag_mass_matrix_vector_relevant = " << diag_mass_matrix_vector_relevant (index_z) -// << ", x = " << point(0) -// << ", y = " << point(1) -// << std::endl; - } - } - } - } - } - distributed_solution.compress(VectorOperation::insert); + if (component == 2) + { + unsigned int index_z = dof_indices[q_point]; + + if (vertex_touched[index_z] == false) + vertex_touched[index_z] = true; + else + continue; + + // the local row where + Point point( + fe_values_face.quadrature_point(q_point)); + + double obstacle_value = obstacle.value(point, 2); + double solution_index_z = solution(index_z); + double gap = obstacle_value - point(2); + + if (lambda(index_z) + / diag_mass_matrix_vector_relevant(index_z) + + c * (solution_index_z - gap) > 0 + && !(constraints_hanging_nodes.is_constrained( + index_z))) + { + constraints.add_line(index_z); + constraints.set_inhomogeneity(index_z, gap); + distributed_solution(index_z) = gap; + + if (locally_owned_dofs.is_element(index_z)) + { + active_set_locally_owned.add_index(index_z); + if (locally_relevant_dofs.is_element(index_z)) + active_set.add_index(index_z); + } + + } + else if (lambda(index_z) + / diag_mass_matrix_vector_relevant(index_z) + + c * (solution_index_z - gap) > 0 + && constraints_hanging_nodes.is_constrained( + index_z)) + { + if (locally_owned_dofs.is_element(index_z)) + { + counter_hanging_nodes += 1; + +// std::cout << "index_z = " << index_z +// << ", lambda = " << lambda (index_z) +// << ", solution_index_z - gap = " << solution_index_z - gap +// << ", diag_mass_matrix_vector_relevant = " << diag_mass_matrix_vector_relevant (index_z) +// << ", x = " << point(0) +// << ", y = " << point(1) +// << std::endl; + } + } + } + } + } + distributed_solution.compress(VectorOperation::insert); - unsigned int sum_contact_constraints = Utilities::MPI::sum( - active_set_locally_owned.n_elements(), mpi_communicator); - pcout << " Size of active set: " << sum_contact_constraints + unsigned int sum_contact_constraints = Utilities::MPI::sum( + active_set_locally_owned.n_elements(), mpi_communicator); + pcout << " Size of active set: " << sum_contact_constraints << std::endl; - unsigned int sum_contact_hanging_nodes = Utilities::MPI::sum( - counter_hanging_nodes, mpi_communicator); - pcout << " Number of hanging nodes in contact: " + unsigned int sum_contact_hanging_nodes = Utilities::MPI::sum( + counter_hanging_nodes, mpi_communicator); + pcout << " Number of hanging nodes in contact: " << sum_contact_hanging_nodes << std::endl; - solution = distributed_solution; + solution = distributed_solution; - constraints.close(); + constraints.close(); - // constraints_dirichlet_hanging_nodes.print (std::cout); + // constraints_dirichlet_hanging_nodes.print (std::cout); - constraints.merge(constraints_dirichlet_hanging_nodes); + constraints.merge(constraints_dirichlet_hanging_nodes); - //constraints.print (std::cout); - } + //constraints.print (std::cout); + } // @sect4{PlasticityContactProblem::dirichlet_constraints} @@ -1459,38 +1459,38 @@ namespace Step42 // the Dirichlet boundary values as well as the // hanging nodes constraints. template - void - PlasticityContactProblem::dirichlet_constraints () - { - /* boundary_indicators: - _______ - / 1 /| - /______ / | - 8| | 8| - | 8 | / - |_______|/ - 6 - */ - - constraints_dirichlet_hanging_nodes.reinit(locally_relevant_dofs); - constraints_dirichlet_hanging_nodes.merge(constraints_hanging_nodes); - - // interpolate all components of the solution - VectorTools::interpolate_boundary_values(dof_handler, - base_mesh == "box" ? 6 : 0, EquationData::BoundaryValues(), - constraints_dirichlet_hanging_nodes, ComponentMask()); - - // interpolate x- and y-components of the - // solution (this is a bit mask, so apply - // operator| ) - FEValuesExtractors::Scalar x_displacement(0); - FEValuesExtractors::Scalar y_displacement(1); - VectorTools::interpolate_boundary_values(dof_handler, 8, - EquationData::BoundaryValues(), - constraints_dirichlet_hanging_nodes, - (fe.component_mask(x_displacement) | fe.component_mask(y_displacement))); - constraints_dirichlet_hanging_nodes.close(); - } + void + PlasticityContactProblem::dirichlet_constraints () + { + /* boundary_indicators: + _______ + / 1 /| + /______ / | + 8| | 8| + | 8 | / + |_______|/ + 6 + */ + + constraints_dirichlet_hanging_nodes.reinit(locally_relevant_dofs); + constraints_dirichlet_hanging_nodes.merge(constraints_hanging_nodes); + + // interpolate all components of the solution + VectorTools::interpolate_boundary_values(dof_handler, + base_mesh == "box" ? 6 : 0, EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, ComponentMask()); + + // interpolate x- and y-components of the + // solution (this is a bit mask, so apply + // operator| ) + FEValuesExtractors::Scalar x_displacement(0); + FEValuesExtractors::Scalar y_displacement(1); + VectorTools::interpolate_boundary_values(dof_handler, 8, + EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, + (fe.component_mask(x_displacement) | fe.component_mask(y_displacement))); + constraints_dirichlet_hanging_nodes.close(); + } // @sect4{PlasticityContactProblem::solve} @@ -1519,60 +1519,60 @@ namespace Step42 // value gamma the linear system becomes // almost semi definite but still symmetric. template - void - PlasticityContactProblem::solve () - { - TimerOutput::Scope t(computing_timer, "Solve"); + void + PlasticityContactProblem::solve () + { + TimerOutput::Scope t(computing_timer, "Solve"); - TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton); - distributed_solution = solution; + TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton); + distributed_solution = solution; - constraints_hanging_nodes.set_zero(distributed_solution); - constraints_hanging_nodes.set_zero(system_rhs_newton); - distributed_solution.compress(VectorOperation::insert); - system_rhs_newton.compress(VectorOperation::insert); + constraints_hanging_nodes.set_zero(distributed_solution); + constraints_hanging_nodes.set_zero(system_rhs_newton); + distributed_solution.compress(VectorOperation::insert); + system_rhs_newton.compress(VectorOperation::insert); - { - TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner"); - preconditioner_u.initialize(system_matrix_newton, additional_data); - } + { + TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner"); + preconditioner_u.initialize(system_matrix_newton, additional_data); + } - { - TimerOutput::Scope t(computing_timer, "Solve: iterate"); - - PrimitiveVectorMemory mem; - TrilinosWrappers::MPI::Vector tmp(system_rhs_newton); - // 1e-4 seems to be the fasted option altogether, but to get more - // reproducible parallel benchmark results, we use a small residual: - double relative_accuracy = 1e-8; - if (output_dir.compare("its/") == 0) - relative_accuracy = 1e-4; - - const double solver_tolerance = relative_accuracy - * system_matrix_newton.residual(tmp, distributed_solution, - system_rhs_newton); - - SolverControl solver_control(system_matrix_newton.m(), - solver_tolerance); - SolverBicgstab solver(solver_control, - mem/*, + { + TimerOutput::Scope t(computing_timer, "Solve: iterate"); + + PrimitiveVectorMemory mem; + TrilinosWrappers::MPI::Vector tmp(system_rhs_newton); + // 1e-4 seems to be the fasted option altogether, but to get more + // reproducible parallel benchmark results, we use a small residual: + double relative_accuracy = 1e-8; + if (output_dir.compare("its/") == 0) + relative_accuracy = 1e-4; + + const double solver_tolerance = relative_accuracy + * system_matrix_newton.residual(tmp, distributed_solution, + system_rhs_newton); + + SolverControl solver_control(system_matrix_newton.m(), + solver_tolerance); + SolverBicgstab solver(solver_control, + mem/*, SolverFGMRES:: AdditionalData(30, true)*/); - solver.solve(system_matrix_newton, distributed_solution, - system_rhs_newton, preconditioner_u); + solver.solve(system_matrix_newton, distributed_solution, + system_rhs_newton, preconditioner_u); - pcout << " Error: " << solver_control.initial_value() - << " -> " << solver_control.last_value() << " in " - << solver_control.last_step() << " Bicgstab iterations." - << std::endl; + pcout << " Error: " << solver_control.initial_value() + << " -> " << solver_control.last_value() << " in " + << solver_control.last_step() << " Bicgstab iterations." + << std::endl; - number_iterations += solver_control.last_step(); - } + number_iterations += solver_control.last_step(); + } - constraints.distribute(distributed_solution); + constraints.distribute(distributed_solution); - solution = distributed_solution; - } + solution = distributed_solution; + } // @sect4{PlasticityContactProblem::solve_newton} @@ -1582,295 +1582,295 @@ namespace Step42 // will be used only if necessary. To obtain a good and reasonable // starting value we solve an elastic problem in very first step (j=1). template - void - PlasticityContactProblem::solve_newton () - { - TimerOutput::Scope t(computing_timer, "solve newton setup"); + void + PlasticityContactProblem::solve_newton () + { + TimerOutput::Scope t(computing_timer, "solve newton setup"); - double resid = 0; - double resid_old = 100000; - TrilinosWrappers::MPI::Vector old_solution(system_rhs_newton); - TrilinosWrappers::MPI::Vector res(system_rhs_newton); - TrilinosWrappers::MPI::Vector tmp_vector(system_rhs_newton); + double resid = 0; + double resid_old = 100000; + TrilinosWrappers::MPI::Vector old_solution(system_rhs_newton); + TrilinosWrappers::MPI::Vector res(system_rhs_newton); + TrilinosWrappers::MPI::Vector tmp_vector(system_rhs_newton); - std::vector < std::vector > constant_modes; - DoFTools::extract_constant_modes(dof_handler, ComponentMask(), - constant_modes); + std::vector < std::vector > constant_modes; + DoFTools::extract_constant_modes(dof_handler, ComponentMask(), + constant_modes); - double sigma_hlp = sigma_0; + double sigma_hlp = sigma_0; - additional_data.constant_modes = constant_modes; - additional_data.elliptic = true; - additional_data.n_cycles = 1; - additional_data.w_cycle = false; - additional_data.output_details = false; - additional_data.smoother_sweeps = 2; - additional_data.aggregation_threshold = 1e-2; + additional_data.constant_modes = constant_modes; + additional_data.elliptic = true; + additional_data.n_cycles = 1; + additional_data.w_cycle = false; + additional_data.output_details = false; + additional_data.smoother_sweeps = 2; + additional_data.aggregation_threshold = 1e-2; - IndexSet active_set_old(active_set); + IndexSet active_set_old(active_set); - t.stop(); // stop newton setup timer + t.stop(); // stop newton setup timer - unsigned int j = 1; - unsigned int number_assemble_system = 0; - for (; j <= 100; j++) - { - if (transfer_solution) - { - if (transfer_solution && j == 1 && cycle == 0) - plast_lin_hard->set_sigma_0(1e+10); - else if (transfer_solution && (j == 2 || cycle > 0)) - plast_lin_hard->set_sigma_0(sigma_hlp); - } - else - { - if (j == 1) - plast_lin_hard->set_sigma_0(1e+10); - else - plast_lin_hard->set_sigma_0(sigma_hlp); - } + unsigned int j = 1; + unsigned int number_assemble_system = 0; + for (; j <= 100; j++) + { + if (transfer_solution) + { + if (transfer_solution && j == 1 && cycle == 0) + plast_lin_hard->set_sigma_0(1e+10); + else if (transfer_solution && (j == 2 || cycle > 0)) + plast_lin_hard->set_sigma_0(sigma_hlp); + } + else + { + if (j == 1) + plast_lin_hard->set_sigma_0(1e+10); + else + plast_lin_hard->set_sigma_0(sigma_hlp); + } - pcout << " " << std::endl; - pcout << " Newton iteration " << j << std::endl; - pcout << " Updating active set..." << std::endl; + pcout << " " << std::endl; + pcout << " Newton iteration " << j << std::endl; + pcout << " Updating active set..." << std::endl; - { - TimerOutput::Scope t(computing_timer, "update active set"); - update_solution_and_constraints(); - } + { + TimerOutput::Scope t(computing_timer, "update active set"); + update_solution_and_constraints(); + } - pcout << " Assembling system... " << std::endl; - system_matrix_newton = 0; - system_rhs_newton = 0; - assemble_nl_system(solution); //compute Newton-Matrix - - number_assemble_system += 1; - - pcout << " Solving system... " << std::endl; - solve(); - - TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton); - distributed_solution = solution; - - // We handle a highly nonlinear problem so we have to damp - // the Newtons method. We refer that we iterate the new solution - // in each Newton step and not only the solution update. - // Since the solution set is a convex set and not a space we - // compute for the damping a linear combination of the - // previous and the current solution to guarantee that the - // damped solution is in our solution set again. - // At most we apply 10 damping steps. - bool damped = false; - tmp_vector = old_solution; - double a = 0; - for (unsigned int i = 0; (i < 5) && (!damped); i++) - { - a = std::pow(0.5, static_cast(i)); - old_solution = tmp_vector; - old_solution.sadd(1 - a, a, distributed_solution); - old_solution.compress(VectorOperation::add); + pcout << " Assembling system... " << std::endl; + system_matrix_newton = 0; + system_rhs_newton = 0; + assemble_nl_system(solution); //compute Newton-Matrix + + number_assemble_system += 1; + + pcout << " Solving system... " << std::endl; + solve(); + + TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton); + distributed_solution = solution; + + // We handle a highly nonlinear problem so we have to damp + // the Newtons method. We refer that we iterate the new solution + // in each Newton step and not only the solution update. + // Since the solution set is a convex set and not a space we + // compute for the damping a linear combination of the + // previous and the current solution to guarantee that the + // damped solution is in our solution set again. + // At most we apply 10 damping steps. + bool damped = false; + tmp_vector = old_solution; + double a = 0; + for (unsigned int i = 0; (i < 5) && (!damped); i++) + { + a = std::pow(0.5, static_cast(i)); + old_solution = tmp_vector; + old_solution.sadd(1 - a, a, distributed_solution); + old_solution.compress(VectorOperation::add); - TimerOutput::Scope t(computing_timer, "Residual and lambda"); + TimerOutput::Scope t(computing_timer, "Residual and lambda"); - system_rhs_newton = 0; - system_rhs_lambda = 0; + system_rhs_newton = 0; + system_rhs_lambda = 0; - solution = old_solution; - residual_nl_system(solution); - res = system_rhs_newton; + solution = old_solution; + residual_nl_system(solution); + res = system_rhs_newton; - const unsigned int start_res = (res.local_range().first), - end_res = (res.local_range().second); - for (unsigned int n = start_res; n < end_res; ++n) - if (constraints.is_inhomogeneously_constrained(n)) - res(n) = 0; + const unsigned int start_res = (res.local_range().first), + end_res = (res.local_range().second); + for (unsigned int n = start_res; n < end_res; ++n) + if (constraints.is_inhomogeneously_constrained(n)) + res(n) = 0; - res.compress(VectorOperation::insert); + res.compress(VectorOperation::insert); - resid = res.l2_norm(); + resid = res.l2_norm(); - if (resid < resid_old) - damped = true; + if (resid < resid_old) + damped = true; - pcout << " Residual of the non-contact part of the system: " + pcout << " Residual of the non-contact part of the system: " << resid << std::endl << " with a damping parameter alpha = " << a << std::endl; - // The previous iteration of step 0 is the solution of an elastic problem. - // So a linear combination of a plastic and an elastic solution makes no sense - // since the elastic solution is not in the convex set of the plastic solution. - if (!transfer_solution && j == 2) - break; - if (transfer_solution && j == 2 && cycle == 0) - break; - } + // The previous iteration of step 0 is the solution of an elastic problem. + // So a linear combination of a plastic and an elastic solution makes no sense + // since the elastic solution is not in the convex set of the plastic solution. + if (!transfer_solution && j == 2) + break; + if (transfer_solution && j == 2 && cycle == 0) + break; + } - resid_old = resid; + resid_old = resid; - resid_vector = system_rhs_lambda; - resid_vector.compress(VectorOperation::insert); + resid_vector = system_rhs_lambda; + resid_vector.compress(VectorOperation::insert); - int is_my_set_changed = (active_set == active_set_old) ? 0 : 1; - int num_changed = Utilities::MPI::sum(is_my_set_changed, - MPI_COMM_WORLD); - if (num_changed == 0) - { - pcout << " Active set did not change!" << std::endl; - if (output_dir.compare("its/") != 0 && resid < 1e-7) - break; - else if (output_dir.compare("its/") == 0 && resid < 1e-10) - break; - } - active_set_old = active_set; - } + int is_my_set_changed = (active_set == active_set_old) ? 0 : 1; + int num_changed = Utilities::MPI::sum(is_my_set_changed, + MPI_COMM_WORLD); + if (num_changed == 0) + { + pcout << " Active set did not change!" << std::endl; + if (output_dir.compare("its/") != 0 && resid < 1e-7) + break; + else if (output_dir.compare("its/") == 0 && resid < 1e-10) + break; + } + active_set_old = active_set; + } - pcout << "" << std::endl << " Number of assembled systems = " + pcout << "" << std::endl << " Number of assembled systems = " << number_assemble_system << std::endl << " Number of Solver-Iterations = " << number_iterations << std::endl; - } + } // @sect3{The refine_grid function} template - void - PlasticityContactProblem::refine_grid () - { - if (refinement_strategy == RefinementStrategy::refine_global) - { - triangulation.refine_global(1); - } - else - { - Vector estimated_error_per_cell( - triangulation.n_active_cells()); - KellyErrorEstimator::estimate(dof_handler, - QGauss(fe.degree + 2), typename FunctionMap::type(), - solution, estimated_error_per_cell); + void + PlasticityContactProblem::refine_grid () + { + if (refinement_strategy == RefinementStrategy::refine_global) + { + triangulation.refine_global(1); + } + else + { + Vector estimated_error_per_cell( + triangulation.n_active_cells()); + KellyErrorEstimator::estimate(dof_handler, + QGauss(fe.degree + 2), typename FunctionMap::type(), + solution, estimated_error_per_cell); - parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( - triangulation, estimated_error_per_cell, 0.3, 0.03); + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( + triangulation, estimated_error_per_cell, 0.3, 0.03); - triangulation.prepare_coarsening_and_refinement(); - if (transfer_solution) - soltrans->prepare_for_coarsening_and_refinement(solution); + triangulation.prepare_coarsening_and_refinement(); + if (transfer_solution) + soltrans->prepare_for_coarsening_and_refinement(solution); - triangulation.execute_coarsening_and_refinement(); - } - } + triangulation.execute_coarsening_and_refinement(); + } + } // @sect3{The move_mesh function} template - void - PlasticityContactProblem::move_mesh ( - const TrilinosWrappers::MPI::Vector &_complete_displacement) const - { - std::vector vertex_touched(triangulation.n_vertices(), false); + void + PlasticityContactProblem::move_mesh ( + const TrilinosWrappers::MPI::Vector &_complete_displacement) const + { + std::vector vertex_touched(triangulation.n_vertices(), false); - for (typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(); cell != dof_handler.end(); ++cell) - if (cell->is_locally_owned()) - for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; - ++v) - { - if (vertex_touched[cell->vertex_index(v)] == false) - { - vertex_touched[cell->vertex_index(v)] = true; + for (typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned()) + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; + ++v) + { + if (vertex_touched[cell->vertex_index(v)] == false) + { + vertex_touched[cell->vertex_index(v)] = true; - Point vertex_displacement; - for (unsigned int d = 0; d < dim; ++d) - { - if (_complete_displacement(cell->vertex_dof_index(v, d)) - != 0) - vertex_displacement[d] = _complete_displacement( - cell->vertex_dof_index(v, d)); - } + Point vertex_displacement; + for (unsigned int d = 0; d < dim; ++d) + { + if (_complete_displacement(cell->vertex_dof_index(v, d)) + != 0) + vertex_displacement[d] = _complete_displacement( + cell->vertex_dof_index(v, d)); + } - cell->vertex(v) += vertex_displacement; - } - } - } + cell->vertex(v) += vertex_displacement; + } + } + } // @sect4{PlasticityContactProblem::output_results} template - void - PlasticityContactProblem::output_results ( - const std::string &title) - { - move_mesh(solution); - - // Calculation of the contact forces - TrilinosWrappers::MPI::Vector lambda(solution); - TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton); - const unsigned int start_res = (resid_vector.local_range().first), - end_res = (resid_vector.local_range().second); - for (unsigned int n = start_res; n < end_res; ++n) - if (constraints.is_inhomogeneously_constrained(n)) - distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n); - distributed_lambda.compress(VectorOperation::insert); - constraints_hanging_nodes.distribute(distributed_lambda); - lambda = distributed_lambda; - TrilinosWrappers::MPI::Vector resid_vector_relevant(solution); - TrilinosWrappers::MPI::Vector distributed_resid_vector(resid_vector); - constraints_hanging_nodes.distribute(distributed_resid_vector); - resid_vector_relevant = distributed_resid_vector; - - DataOut data_out; - - data_out.attach_dof_handler(dof_handler); - - const std::vector data_component_interpretation( - dim, DataComponentInterpretation::component_is_part_of_vector); - data_out.add_data_vector(solution, - std::vector < std::string > (dim, "Displacement"), - DataOut::type_dof_data, data_component_interpretation); - data_out.add_data_vector(lambda, - std::vector < std::string > (dim, "ContactForce"), - DataOut::type_dof_data, data_component_interpretation); - data_out.add_data_vector(active_set, - std::vector < std::string > (dim, "ActiveSet"), - DataOut::type_dof_data, data_component_interpretation); - data_out.add_data_vector(resid_vector_relevant, - std::vector < std::string > (dim, "Residual"), - DataOut::type_dof_data, data_component_interpretation); - - Vector subdomain(triangulation.n_active_cells()); - for (unsigned int i = 0; i < subdomain.size(); ++i) - subdomain(i) = triangulation.locally_owned_subdomain(); - data_out.add_data_vector(subdomain, "subdomain"); - - data_out.add_data_vector(cell_constitution, "CellConstitution"); - - data_out.build_patches(); - - const std::string filename = - (output_dir + title + "-" - + 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 + title << ".pvtu" << std::endl; - - if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) - { - std::vector < std::string > filenames; - for (unsigned int i = 0; - i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i) - filenames.push_back( - title + "-" + Utilities::int_to_string(i, 4) + ".vtu"); - - std::ofstream master_output((output_dir + title + ".pvtu").c_str()); - data_out.write_pvtu_record(master_output, filenames); - } + void + PlasticityContactProblem::output_results ( + const std::string &title) + { + move_mesh(solution); + + // Calculation of the contact forces + TrilinosWrappers::MPI::Vector lambda(solution); + TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton); + const unsigned int start_res = (resid_vector.local_range().first), + end_res = (resid_vector.local_range().second); + for (unsigned int n = start_res; n < end_res; ++n) + if (constraints.is_inhomogeneously_constrained(n)) + distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n); + distributed_lambda.compress(VectorOperation::insert); + constraints_hanging_nodes.distribute(distributed_lambda); + lambda = distributed_lambda; + TrilinosWrappers::MPI::Vector resid_vector_relevant(solution); + TrilinosWrappers::MPI::Vector distributed_resid_vector(resid_vector); + constraints_hanging_nodes.distribute(distributed_resid_vector); + resid_vector_relevant = distributed_resid_vector; + + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + + const std::vector data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_out.add_data_vector(solution, + std::vector < std::string > (dim, "Displacement"), + DataOut::type_dof_data, data_component_interpretation); + data_out.add_data_vector(lambda, + std::vector < std::string > (dim, "ContactForce"), + DataOut::type_dof_data, data_component_interpretation); + data_out.add_data_vector(active_set, + std::vector < std::string > (dim, "ActiveSet"), + DataOut::type_dof_data, data_component_interpretation); + data_out.add_data_vector(resid_vector_relevant, + std::vector < std::string > (dim, "Residual"), + DataOut::type_dof_data, data_component_interpretation); + + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = triangulation.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + data_out.add_data_vector(cell_constitution, "CellConstitution"); + + data_out.build_patches(); + + const std::string filename = + (output_dir + title + "-" + + 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 + title << ".pvtu" << std::endl; + + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + std::vector < std::string > filenames; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i) + filenames.push_back( + title + "-" + Utilities::int_to_string(i, 4) + ".vtu"); + + std::ofstream master_output((output_dir + title + ".pvtu").c_str()); + data_out.write_pvtu_record(master_output, filenames); + } - TrilinosWrappers::MPI::Vector tmp(solution); - tmp *= -1; - move_mesh(tmp); - } + TrilinosWrappers::MPI::Vector tmp(solution); + tmp *= -1; + move_mesh(tmp); + } // @sect4{PlasticityContactProblem::output_contact_force} @@ -1891,192 +1891,192 @@ namespace Step42 // a processor does not own the cell with the point we have to // catch these cases. template - void - PlasticityContactProblem::output_contact_force ( - const unsigned int cycle) + void + PlasticityContactProblem::output_contact_force ( + const unsigned int cycle) + { + Functions::FEFieldFunction, + TrilinosWrappers::MPI::Vector> solution_function(dof_handler, + solution); + std::cout.precision(10); + + Vector solution_p1(dim); + std::vector > solution_gradient_p1(dim); + + // Here we calculate the contact pressure as a vector lambda. + // If a dof is element of the active set lambda contains the + // nonlinear residual this dof divided by the according entry + // of the mass matrix. In all other dofs lambda will be set to + // zero. + TrilinosWrappers::MPI::Vector lambda(solution); + TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton); + const unsigned int start_res = (resid_vector.local_range().first), + end_res = (resid_vector.local_range().second); + for (unsigned int n = start_res; n < end_res; ++n) + if (constraints.is_inhomogeneously_constrained(n)) + distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n); + else + distributed_lambda(n) = 0; + distributed_lambda.compress(VectorOperation::insert); + constraints_hanging_nodes.distribute(distributed_lambda); + lambda = distributed_lambda; + Functions::FEFieldFunction, + TrilinosWrappers::MPI::Vector> lambda_function(dof_handler, lambda); + + // Here we try to find the MPI-process which owns the cell + // with the point_of_interest. If it is the wrong MPI-process + // we catch this case and set point_found to false. + const Point point_of_interest(0.49, 0.5001, 1.0); + Vector contact_pressure_in_point(dim); + bool point_found = true; + + MPI_Barrier(MPI_COMM_WORLD); + try + { + lambda_function.vector_value(point_of_interest, + contact_pressure_in_point); + } + catch (const typename Functions::FEFieldFunction, + TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &) + { + point_found = false; + } + + if (point_found == true) + { + std::cout << "PoI contact pressure: " << contact_pressure_in_point(2) + << std::endl; + } + + // To obtain the contact force we have to compute an integral of the contact pressure + // in z-direction over the whole contact area. To be accurate enough we use the + // Gaussian quadrature rule with fe.degree + 1. + double contact_force = 0.0; { - Functions::FEFieldFunction, - TrilinosWrappers::MPI::Vector> solution_function(dof_handler, - solution); - std::cout.precision(10); - - Vector solution_p1(dim); - std::vector > solution_gradient_p1(dim); - - // Here we calculate the contact pressure as a vector lambda. - // If a dof is element of the active set lambda contains the - // nonlinear residual this dof divided by the according entry - // of the mass matrix. In all other dofs lambda will be set to - // zero. - TrilinosWrappers::MPI::Vector lambda(solution); - TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton); - const unsigned int start_res = (resid_vector.local_range().first), - end_res = (resid_vector.local_range().second); - for (unsigned int n = start_res; n < end_res; ++n) - if (constraints.is_inhomogeneously_constrained(n)) - distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n); - else - distributed_lambda(n) = 0; - distributed_lambda.compress(VectorOperation::insert); - constraints_hanging_nodes.distribute(distributed_lambda); - lambda = distributed_lambda; - Functions::FEFieldFunction, - TrilinosWrappers::MPI::Vector> lambda_function(dof_handler, lambda); - - // Here we try to find the MPI-process which owns the cell - // with the point_of_interest. If it is the wrong MPI-process - // we catch this case and set point_found to false. - const Point point_of_interest(0.49, 0.5001, 1.0); - Vector contact_pressure_in_point(dim); - bool point_found = true; - - MPI_Barrier(MPI_COMM_WORLD); - try - { - lambda_function.vector_value(point_of_interest, - contact_pressure_in_point); - } - catch (const typename Functions::FEFieldFunction, - TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &) - { - point_found = false; - } + QGauss face_quadrature_formula(fe.degree + 1); - if (point_found == true) - { - std::cout << "PoI contact pressure: " << contact_pressure_in_point(2) - << std::endl; - } + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); - // To obtain the contact force we have to compute an integral of the contact pressure - // in z-direction over the whole contact area. To be accurate enough we use the - // Gaussian quadrature rule with fe.degree + 1. - double contact_force = 0.0; - { - QGauss face_quadrature_formula(fe.degree + 1); + const unsigned int n_face_q_points = face_quadrature_formula.size(); - FEFaceValues fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points | update_JxW_values); + const FEValuesExtractors::Vector displacement(0); - const unsigned int n_face_q_points = face_quadrature_formula.size(); + 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 = 0; + face < GeometryInfo::faces_per_cell; ++face) + if (cell->face(face)->at_boundary() + && cell->face(face)->boundary_indicator() == 1) + { + fe_values_face.reinit(cell, face); - const FEValuesExtractors::Vector displacement(0); + std::vector > lambda_values(n_face_q_points); + fe_values_face[displacement].get_function_values(lambda, + lambda_values); - 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 = 0; - face < GeometryInfo::faces_per_cell; ++face) - if (cell->face(face)->at_boundary() - && cell->face(face)->boundary_indicator() == 1) + for (unsigned int q_point = 0; q_point < n_face_q_points; + ++q_point) { - fe_values_face.reinit(cell, face); - - std::vector > lambda_values(n_face_q_points); - fe_values_face[displacement].get_function_values(lambda, - lambda_values); - - for (unsigned int q_point = 0; q_point < n_face_q_points; - ++q_point) - { - contact_force += lambda_values[q_point][2] - * fe_values_face.JxW(q_point); - } + contact_force += lambda_values[q_point][2] + * fe_values_face.JxW(q_point); } - contact_force = Utilities::MPI::sum(contact_force, MPI_COMM_WORLD); - pcout << "Contact force = " << contact_force << std::endl; - } - MPI_Barrier(MPI_COMM_WORLD); + } + contact_force = Utilities::MPI::sum(contact_force, MPI_COMM_WORLD); + pcout << "Contact force = " << contact_force << std::endl; } + MPI_Barrier(MPI_COMM_WORLD); + } // @sect4{PlasticityContactProblem::run} template - void - PlasticityContactProblem::run () - { + void + PlasticityContactProblem::run () + { - if (obstacle_filename != "") - { - pcout << "Read the obstacle from '" << obstacle_filename << "' ... " + if (obstacle_filename != "") + { + pcout << "Read the obstacle from '" << obstacle_filename << "' ... " << std::flush; - input_obstacle.reset(new Input(obstacle_filename.c_str())); - pcout << "done." << std::endl; - } + input_obstacle.reset(new Input(obstacle_filename.c_str())); + pcout << "done." << std::endl; + } - computing_timer.reset(); - for (cycle = 0; cycle < n_cycles; ++cycle) + computing_timer.reset(); + for (cycle = 0; cycle < n_cycles; ++cycle) + { { - { - TimerOutput::Scope t(computing_timer, "Setup"); + TimerOutput::Scope t(computing_timer, "Setup"); - pcout << std::endl; - pcout << "Cycle " << cycle << ':' << std::endl; - - if (cycle == 0) - { - make_grid(); - } - else - { - TimerOutput::Scope t(computing_timer, "Setup: refine mesh"); - if (transfer_solution) - soltrans.reset( - new parallel::distributed::SolutionTransfer(dof_handler)); - refine_grid(); - } + pcout << std::endl; + pcout << "Cycle " << cycle << ':' << std::endl; - setup_system(); + if (cycle == 0) + { + make_grid(); + } + else + { + TimerOutput::Scope t(computing_timer, "Setup: refine mesh"); + if (transfer_solution) + soltrans.reset( + new parallel::distributed::SolutionTransfer(dof_handler)); + refine_grid(); + } - if (transfer_solution && cycle > 0) - { - TrilinosWrappers::MPI::Vector distributed_solution( - system_rhs_newton); - distributed_solution = solution; - soltrans->interpolate(distributed_solution); - solution = distributed_solution; - residual_nl_system(solution); - resid_vector = system_rhs_lambda; - resid_vector.compress(VectorOperation::insert); - } + setup_system(); + if (transfer_solution && cycle > 0) + { + TrilinosWrappers::MPI::Vector distributed_solution( + system_rhs_newton); + distributed_solution = solution; + soltrans->interpolate(distributed_solution); + solution = distributed_solution; + residual_nl_system(solution); + resid_vector = system_rhs_lambda; + resid_vector.compress(VectorOperation::insert); } - solve_newton(); + } - if (true) //Utilities::MPI::n_mpi_processes(mpi_communicator) <= 64) - { - pcout << " Writing graphical output... " << std::flush; + solve_newton(); - TimerOutput::Scope t(computing_timer, "Graphical output"); + if (true) //Utilities::MPI::n_mpi_processes(mpi_communicator) <= 64) + { + pcout << " Writing graphical output... " << std::flush; - std::ostringstream filename_solution; - filename_solution << "solution-"; - filename_solution << Utilities::int_to_string(cycle, 2); - output_results(filename_solution.str()); - } + TimerOutput::Scope t(computing_timer, "Graphical output"); + + std::ostringstream filename_solution; + filename_solution << "solution-"; + filename_solution << Utilities::int_to_string(cycle, 2); + output_results(filename_solution.str()); + } - computing_timer.print_summary(); - computing_timer.reset(); + computing_timer.print_summary(); + computing_timer.reset(); - Utilities::System::MemoryStats stats; - Utilities::System::get_memory_stats(stats); - pcout << "VMPEAK, Resident in kB: " << stats.VmSize << " " + Utilities::System::MemoryStats stats; + Utilities::System::get_memory_stats(stats); + pcout << "VMPEAK, Resident in kB: " << stats.VmSize << " " << stats.VmRSS << std::endl; - if (base_mesh == "box") - output_contact_force(cycle); - } - } + if (base_mesh == "box") + output_contact_force(cycle); + } + } } // @sect3{The main function} int main ( - int argc, char *argv[]) + int argc, char *argv[]) { using namespace dealii; using namespace Step42; @@ -2092,10 +2092,10 @@ main ( prm.read_input(argv[1]); Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); - { - PlasticityContactProblem<3> problem(prm); - problem.run(); - } + { + PlasticityContactProblem<3> problem(prm); + problem.run(); + } return 0; } -- 2.39.5