From 0031ff1c0177c7fed973ae2620291ad7b7049b6b Mon Sep 17 00:00:00 2001 From: frohne Date: Tue, 27 Aug 2013 12:02:59 +0000 Subject: [PATCH] more comments git-svn-id: https://svn.dealii.org/trunk@30500 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 3498 ++++++++++++++------------- 1 file changed, 1757 insertions(+), 1741 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index ada66613fb..e34ec1fca3 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -21,7 +21,6 @@ * Timo Heister, Texas A&M University, 2013 */ - // @sect3{Include files} // We are using the the same // include files as in step-41: @@ -80,8 +79,9 @@ #include -namespace Step42 { -using namespace dealii; +namespace Step42 +{ + using namespace dealii; // @sect3{The Input class template} @@ -110,135 +110,158 @@ using namespace dealii; // update_solution_and_constraints () of // the class PlasticityContactProblem. -template -class Input { -public: - Input(const char* _name) : - name(_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); - } - - double - hv(int i, int j); - - double - obstacle_function(double x, double y); - - void - read_obstacle(const char* name); - -private: - const char* name; - MPI_Comm mpi_communicator; - ConditionalOStream pcout; - std::vector obstacle_data; - double hx, hy; - int nx, ny; -}; + template + class Input + { + public: + Input ( + const char* _name) + : + name(_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); + } + + double + hv ( + int i, int j); + + double + obstacle_function ( + double x, double y); + + void + read_obstacle ( + const char* name); + + private: + const char* name; + 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(int i, 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 -} + template + double + Input::hv ( + int i, 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(double x, double y) { - int ix = (int) (x / hx); - int iy = (int) (y / hy); - - if (ix < 0) - ix = 0; - - if (iy < 0) - iy = 0; - - if (ix >= nx - 1) - ix = nx - 2; - - if (iy >= ny - 1) - iy = ny - 2; - - double val = 0.0; - { - FullMatrix H(4, 4); - Vector X(4); - Vector b(4); - - double xx = 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; -} + template + double + Input::obstacle_function ( + double x, double y) + { + int ix = (int) (x / hx); + int iy = (int) (y / hy); + + if (ix < 0) + ix = 0; + + if (iy < 0) + iy = 0; + + if (ix >= nx - 1) + ix = nx - 2; + + if (iy >= ny - 1) + iy = ny - 2; + + double val = 0.0; + { + FullMatrix H(4, 4); + Vector X(4); + Vector b(4); + + double xx = 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 char* name) { - std::ifstream f(name); + template + void + Input::read_obstacle ( + const char* name) + { + std::ifstream f(name); - 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 " << ny - << std::endl; -} + pcout << "Resolution of the scanned obstacle picture: " << nx << " x " + << ny << std::endl; + } // @sect3{The ConstitutiveLaw class template} @@ -249,41 +272,48 @@ void Input::read_obstacle(const char* name) { // isotropic hardening. // For gamma = 0 we obtain perfect elastoplastic // behavior. -template -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) { - 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; -}; + template + 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) + { + 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; + }; // The constructor of the ConstitutiveLaw class sets the // required material parameter for our deformable body: @@ -294,37 +324,46 @@ private: // Also it supplies the stress strain tensor of forth order // 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); -} + 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); + } // @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; + 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; - 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} @@ -335,35 +374,40 @@ inline SymmetricTensor<2, dim> ConstitutiveLaw::get_strain( // points. We need this function to calculate the nonlinear // 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; - - SymmetricTensor < 2, dim > deviator_stress_tensor = deviator( - stress_tensor); - - 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; - - stress_strain_tensor += stress_strain_tensor_kappa; - } -} + 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; + + SymmetricTensor<2, dim> deviator_stress_tensor = deviator( + stress_tensor); + + 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; + + stress_strain_tensor += stress_strain_tensor_kappa; + } + } // @sect3{ConstitutiveLaw::linearized_plast_linear_hardening} @@ -377,116 +421,139 @@ void ConstitutiveLaw::plast_linear_hardening( // See // PlasticityContactProblem::assemble_nl_system(TrilinosWrappers::MPI::Vector &u) // where this function is used. -template -void ConstitutiveLaw::linearized_plast_linear_hardening( - SymmetricTensor<4, dim> &stress_strain_tensor_linearized, - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor) { - 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); - - 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_kappa; - stress_strain_tensor_linearized += stress_strain_tensor_kappa; - } -} + 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; + + SymmetricTensor<2, dim> deviator_stress_tensor = deviator( + stress_tensor); + + 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_kappa; + stress_strain_tensor_linearized += stress_strain_tensor_kappa; + } + } -namespace EquationData { + namespace EquationData + { // It possible to apply an additional body force // but in here it is set to zero. -template -class RightHandSide: public Function { -public: - RightHandSide() : - Function(dim) { - } - - virtual double - value(const Point &p, const unsigned int component = 0) 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; - - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - return_value = 0.0; - - return return_value; -} + template + class RightHandSide : public Function + { + public: + RightHandSide () + : + Function(dim) + { + } -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); -} + virtual double + value ( + const Point &p, const unsigned int component = 0) 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; + + if (component == 0) + return_value = 0.0; + if (component == 1) + return_value = 0.0; + if (component == 2) + return_value = 0.0; + + return return_value; + } + + template + 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 { -public: - BoundaryValues() : - Function(dim) { - } - ; - - virtual double - value(const Point &p, const unsigned int component = 0) 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; - - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - return_value = 0.0; - - return return_value; -} - -template -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); -} + template + class BoundaryValues : public Function + { + public: + BoundaryValues () + : + Function(dim) + { + } + ; + + virtual double + value ( + const Point &p, const unsigned int component = 0) 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; + + if (component == 0) + return_value = 0.0; + if (component == 1) + return_value = 0.0; + if (component == 2) + return_value = 0.0; + + return return_value; + } + + template + 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 @@ -494,60 +561,72 @@ void BoundaryValues::vector_value(const Point &p, // 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 -template -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){ - } - - virtual double - value(const Point &p, const unsigned int component = 0) 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; -}; - -template -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; - } -} + template + 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) + { + } -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); -} -} + virtual double + value ( + const Point &p, const unsigned int component = 0) 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; + }; + + template + 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; + } + } + + 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); + } + } // @sect3{The PlasticityContactProblem class template} @@ -566,116 +645,124 @@ void Obstacle::vector_value(const Point &p, // situation and to handle the nonlinear // operator for the constitutive law. -template -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 - { - 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; -}; + template + 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 + { + 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} @@ -683,655 +770,658 @@ private: // template that makes use of the functions // 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 if (strat == "fix dofs") - refinement_strategy = RefinementStrategy::refine_fix_dofs; - 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; -} + 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 " + << (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'"); - -} + 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'"); + + } Point<3> - rotate_half_sphere(const Point<3> &in) + rotate_half_sphere ( + const Point<3> &in) { return Point<3>(in(2), in(1), -in(0)); } // @sect4{PlasticityContactProblem::make_grid} - -template -void PlasticityContactProblem::make_grid() { - if (base_mesh == "half sphere") + template + void + PlasticityContactProblem::make_grid () { - Point < dim > center(0, 0, 0); - double radius = 0.8; - GridGenerator::half_hyper_ball(triangulation, center, radius); - GridTools::transform(&rotate_half_sphere, triangulation); - Point < dim > 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); + 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); + + GridGenerator::hyper_rectangle(triangulation, p1, p2); to_refine_factor = 0.3; to_coarsen_factor = 0.03; - return; + + Triangulation<3>::active_cell_iterator cell = + triangulation.begin_active(), endc = triangulation.end(); + + /* 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); } - - Point < dim > p1(0, 0, 0); - Point < dim > p2(1.0, 1.0, 1.0); - unsigned int ref = n_initial_refinements; - if (refinement_strategy == RefinementStrategy::refine_fix_dofs) + template + void + PlasticityContactProblem::setup_system () { - /** - * This complicated logic creates a mesh and a refinement fraction to_refine_factor, - * so that the resulting mesh after adaptive refinement has approximately - * 2^n_refinements_global*300 dofs. This allows parallel scalability tests. - * About 5%-10% of the cells are being adaptively refined. - * We start with a 3x3,4x4, or 5x5 base mesh (whichever is closed in cell - * count). - */ - unsigned int ref = (n_initial_refinements + 1) / 3; - unsigned int remain = n_initial_refinements + 1 - ref * 3; - unsigned int rep = 3; - if (remain == 1) - rep = 4; - else if (remain == 2) - rep = 5; - - unsigned int n_cells_x = (1 << ref) * rep; - unsigned int goal_dofs = (1 << n_initial_refinements) * 300; - double goal_cells = std::pow(std::pow(goal_dofs / 3.0, 1.0 / 3.0) - 1.0, - 3.0); - double n_cells = std::pow(n_cells_x, 3.0); - to_refine_factor = (goal_cells - n_cells) / n_cells; - //convert from fraction of cells to add to fraction of cells to refine: - to_refine_factor /= 7.0; - to_coarsen_factor = 0.0; - - std::vector repet(3); - repet[0] = rep; - repet[1] = rep; - repet[2] = rep; - - GridGenerator::subdivided_hyper_rectangle(triangulation, repet, p1, p2); + // setup dofs + { + 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); + } + + // 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()); + } + + // 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)); + + sp.compress(); + + 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); + + number_iterations = 0; + + diag_mass_matrix_vector.compress(VectorOperation::insert); + + // remove the mass matrix entries from the matrix: + mass_matrix = 0; + } } - else + + template + void + PlasticityContactProblem::assemble_nl_system ( + TrilinosWrappers::MPI::Vector &u) { - GridGenerator::hyper_rectangle(triangulation, p1, p2); - to_refine_factor = 0.3; - to_coarsen_factor = 0.03; + TimerOutput::Scope t(computing_timer, "Assembling"); + + 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)); + + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const EquationData::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); + + std::vector local_dof_indices(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); + + 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; + + 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); + + 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]); + + 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)); + + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const EquationData::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)); + + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const FEValuesExtractors::Vector displacement(0); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); + + unsigned int elast_points = 0; + unsigned int plast_points = 0; + double yield = 0; + unsigned int cell_number = 0; + cell_constitution = 0; + + for (; cell != endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_rhs = 0; + + 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); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + SymmetricTensor<4, dim> stress_strain_tensor; + SymmetricTensor<2, dim> stress_tensor; + + plast_lin_hard->plast_linear_hardening(stress_strain_tensor, + strain_tensor[q_point], elast_points, plast_points, yield); + + 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->get_dof_indices(local_dof_indices); + constraints_dirichlet_hanging_nodes.distribute_local_to_global( + cell_rhs, local_dof_indices, system_rhs_newton); + + for (unsigned int i = 0; i < dofs_per_cell; i++) + system_rhs_lambda(local_dof_indices[i]) += cell_rhs(i); + + cell_number += 1; + } + else + { + cell_constitution(cell_number) = 0; + cell_number += 1; + }; - Triangulation<3>::active_cell_iterator cell = triangulation.begin_active(), - endc = triangulation.end(); - - /* boundary_indicators: - _______ - / 1 /| - /______ / | - 8| | 8| - | 8 | / - |_______|/ - 6 - */ - - for (; cell != endc; ++cell) - for (unsigned int face = 0; - face < GeometryInfo < dim > ::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(ref); -} - -template -void PlasticityContactProblem::setup_system() { - // setup dofs - { - 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); - } - - // 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()); - } - - // 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)); - - sp.compress(); - - 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); - - number_iterations = 0; - - diag_mass_matrix_vector.compress(VectorOperation::insert); - - // 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"); - - QGauss quadrature_formula(fe.degree + 1); - QGauss face_quadrature_formula(fe.degree + 1); - - FEValues < dim - > fe_values(fe, quadrature_formula, - UpdateFlags( - update_values | update_gradients | update_q_points - | update_JxW_values)); - - FEFaceValues < dim - > fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points - | update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - const EquationData::RightHandSide right_hand_side; - std::vector < Vector - > right_hand_side_values(n_q_points, Vector(dim)); - std::vector < 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); - - std::vector local_dof_indices(dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); - - 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; - - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - right_hand_side_values); - - std::vector < SymmetricTensor<2, dim> > 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; - - 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 face = 0; - face < GeometryInfo < dim > ::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 < dim - > fe_values(fe, quadrature_formula, - UpdateFlags( - update_values | update_gradients | update_q_points - | update_JxW_values)); - - FEFaceValues < dim - > fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points - | update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - const EquationData::RightHandSide right_hand_side; - std::vector < Vector - > right_hand_side_values(n_q_points, Vector(dim)); - std::vector < Vector - > right_hand_side_values_face(n_face_q_points, Vector(dim)); - - Vector cell_rhs(dofs_per_cell); - - std::vector local_dof_indices(dofs_per_cell); - - const FEValuesExtractors::Vector displacement(0); - - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); - - unsigned int elast_points = 0; - unsigned int plast_points = 0; - double yield = 0; - unsigned int cell_number = 0; - cell_constitution = 0; - - for (; cell != endc; ++cell) - if (cell->is_locally_owned()) { - fe_values.reinit(cell); - cell_rhs = 0; - - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - right_hand_side_values); - - std::vector < SymmetricTensor<2, dim> > 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; - SymmetricTensor < 2, dim > stress_tensor; - - plast_lin_hard->plast_linear_hardening(stress_strain_tensor, - strain_tensor[q_point], elast_points, plast_points, - yield); - - 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 < dim > ::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_dirichlet_hanging_nodes.distribute_local_to_global( - cell_rhs, local_dof_indices, system_rhs_newton); - - for (unsigned int i=0; i -void PlasticityContactProblem::assemble_mass_matrix_diagonal( - TrilinosWrappers::SparseMatrix &mass_matrix) { - QGaussLobatto < dim - 1 > face_quadrature_formula(fe.degree + 1); + template + void + PlasticityContactProblem::assemble_mass_matrix_diagonal ( + TrilinosWrappers::SparseMatrix &mass_matrix) + { + QGaussLobatto face_quadrature_formula(fe.degree + 1); - FEFaceValues < dim - > 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 < dim > ::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); - for (unsigned int i=0; i -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 < dim > ::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 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; + 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(); + + // 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; + + 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) @@ -1340,32 +1430,32 @@ void PlasticityContactProblem::update_solution_and_constraints() { // << ", x = " << point(0) // << ", y = " << point(1) // << std::endl; - } - } - } - } - } - distributed_solution.compress(VectorOperation::insert); + } + } + } + } + } + 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 - << std::endl; - 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; + 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: " + << 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} @@ -1373,38 +1463,39 @@ void PlasticityContactProblem::update_solution_and_constraints() { // constraints_dirichlet_hanging_nodes. It contains // 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(); -} + 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(); + } // @sect4{PlasticityContactProblem::solve} @@ -1432,59 +1523,61 @@ void PlasticityContactProblem::dirichlet_constraints() { // instead of CG. For a very small hardening // value gamma the linear system becomes // almost semi definite but still symmetric. -template -void PlasticityContactProblem::solve() { - TimerOutput::Scope t(computing_timer, "Solve"); - - 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); - - { - TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner"); - preconditioner_u.initialize(system_matrix_newton, additional_data); - } - - { - TimerOutput::Scope t(computing_timer, "Solve: iterate"); - - PrimitiveVectorMemory < TrilinosWrappers::MPI::Vector > 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 < TrilinosWrappers::MPI::Vector - > solver(solver_control, mem/*, - SolverFGMRES:: - AdditionalData(30, true)*/); - 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; - - number_iterations += solver_control.last_step(); - } - - constraints.distribute(distributed_solution); - - solution = distributed_solution; -} + template + void + PlasticityContactProblem::solve () + { + TimerOutput::Scope t(computing_timer, "Solve"); + + 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); + + { + 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/*, + SolverFGMRES:: + AdditionalData(30, true)*/); + 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; + + number_iterations += solver_control.last_step(); + } + + constraints.distribute(distributed_solution); + + solution = distributed_solution; + } // @sect4{PlasticityContactProblem::solve_newton} @@ -1493,281 +1586,296 @@ void PlasticityContactProblem::solve() { // iteration and the inner loop for the damping steps which // 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"); - - 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); - - 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; - - IndexSet active_set_old(active_set); - - 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); - } - - 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(); - } - - 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"); - - system_rhs_newton = 0; - system_rhs_lambda = 0; - - 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; - - res.compress(VectorOperation::insert); - - resid = res.l2_norm(); - - if (resid < resid_old) - damped = true; - - 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; - } - - resid_old = resid; - - 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!" < > constant_modes; + DoFTools::extract_constant_modes(dof_handler, ComponentMask(), + constant_modes); + + 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; + + IndexSet active_set_old(active_set); + + 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); + } + + 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(); + } + + 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"); + + system_rhs_newton = 0; + system_rhs_lambda = 0; + + 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; + + res.compress(VectorOperation::insert); + + resid = res.l2_norm(); + + if (resid < resid_old) + damped = true; + + 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; + } + + resid_old = resid; + + 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; + } + + pcout << "" << std::endl << " Number of assembled systems = " + << number_assemble_system << std::endl + << " Number of Solver-Iterations = " << number_iterations + << std::endl; } - else - { - Vector estimated_error_per_cell(triangulation.n_active_cells()); - KellyErrorEstimator < dim - > ::estimate(dof_handler, QGauss < dim - 1 > (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); - - triangulation.prepare_coarsening_and_refinement(); - if (transfer_solution) - soltrans->prepare_for_coarsening_and_refinement(solution); +// @sect3{The refine_grid function} - triangulation.execute_coarsening_and_refinement(); + 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); + + 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.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); - - 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 < dim > ::vertices_per_cell; ++v) { - if (vertex_touched[cell->vertex_index(v)] == false) { - vertex_touched[cell->vertex_index(v)] = true; - - Point < dim > 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; - } - } -} + template + 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; + + 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; + } + } + } // @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 < dim > 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 < dim > ::type_dof_data, data_component_interpretation); - data_out.add_data_vector(lambda, - std::vector < std::string > (dim, "ContactForce"), - DataOut < dim > ::type_dof_data, data_component_interpretation); - data_out.add_data_vector(active_set, - std::vector < std::string > (dim, "ActiveSet"), - DataOut < dim > ::type_dof_data, data_component_interpretation); - data_out.add_data_vector(resid_vector_relevant, - std::vector < std::string > (dim, "Residual"), - DataOut < dim > ::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 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); -} + 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); + } + + TrilinosWrappers::MPI::Vector tmp(solution); + tmp *= -1; + move_mesh(tmp); + } // @sect4{PlasticityContactProblem::output_contact_force} @@ -1782,309 +1890,217 @@ void PlasticityContactProblem::output_results( // it is important to apply contraints_hanging_nodes.distribute // to the distributed_lambda vector. // To calculate the contact pressure in a certain point in the -// contact area, we have make use of the Functions::FEFieldFunction -// In parallel this is little tricky because we have to find the -// process with the right cell which contains this point. -template -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 < Tensor<1, dim> > solution_gradient_p1(dim); - - const Point p1_of_interest(0.5001, 0.5001, 0.9501); - bool point1_found = true; - bool point2_found = true; - - // 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); - 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); - const Point p2_of_interest(0.49, 0.5001, 1.0); - Vector lambda_p2(dim); - - MPI_Barrier(MPI_COMM_WORLD); - try { - lambda_function.vector_value(p2_of_interest, lambda_p2); - } catch (const typename Functions::FEFieldFunction, - TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &) { - point2_found = false; - } - - if (point2_found == true) { - std::cout << "PoI lambda_z: " << lambda_p2(2) << std::endl; - } - - // Integral of the contact force in z-direction over the whole contact area. - double contact_force = 0.0; - { - QGauss< dim - 1 > face_quadrature_formula(fe.degree + 1); - - FEFaceValues < dim > fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points - | update_JxW_values); - - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - const FEValuesExtractors::Vector displacement(0); - - 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 < dim > ::faces_per_cell; ++face) - if (cell->face(face)->at_boundary() - && cell->face(face)->boundary_indicator() == 1) { - fe_values_face.reinit(cell, face); - - std::vector < Tensor<1, dim> - > 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 = Utilities::MPI::sum(contact_force, - MPI_COMM_WORLD); - pcout << "Contact force = " << contact_force << std::endl; - } - - // To calculate the contact area between deformable body and obstacle - double contact_area = 0.0; - { - move_mesh(solution); - - QGaussLobatto < dim - 1 > face_quadrature_formula(fe.degree + 1); - - FEFaceValues < dim > fe_values_face(fe, face_quadrature_formula, - update_values | update_quadrature_points - | update_JxW_values); - - const unsigned int dofs_per_face = fe.dofs_per_face; - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - const FEValuesExtractors::Vector displacement(0); - 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 < dim > ::faces_per_cell; ++face){ - if (cell->face(face)->at_boundary() - && cell->face(face)->boundary_indicator() == 1) { - fe_values_face.reinit(cell, face); - - unsigned int contact_counter = 0; - std::vector dof_indices (dofs_per_face); - cell->face(face)->get_dof_indices (dof_indices); - - for (unsigned int q_point=0; q_point + 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; + { + QGauss face_quadrature_formula(fe.degree + 1); + + FEFaceValues fe_values_face(fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); + + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const FEValuesExtractors::Vector displacement(0); + + 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); + + 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 = 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() { - - 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; - } + template + void + PlasticityContactProblem::run () + { - computing_timer.reset(); - for (cycle = 0; cycle < n_cycles; ++cycle) { - { - 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(); - } - - 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; - - 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(); - - 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_for_benchmark(cycle); - } -} + 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; + } + + computing_timer.reset(); + for (cycle = 0; cycle < n_cycles; ++cycle) + { + { + 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(); + } + + 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; + + 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(); + + 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); + } + } } // @sect3{The main function} -int main(int argc, char *argv[]) { - using namespace dealii; - using namespace Step42; - - deallog.depth_console(0); - ParameterHandler prm; - PlasticityContactProblem<3>::declare(prm); - if (argc!=2) - { - prm.print_parameters(std::cout, ParameterHandler::Text); - return 0; - } - - prm.read_input(argv[1]); - Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); - { - PlasticityContactProblem<3> problem(prm); - problem.run(); - } - - return 0; +int +main ( + int argc, char *argv[]) +{ + using namespace dealii; + using namespace Step42; + + deallog.depth_console(0); + ParameterHandler prm; + PlasticityContactProblem<3>::declare(prm); + if (argc != 2) + { + prm.print_parameters(std::cout, ParameterHandler::Text); + return 0; + } + + prm.read_input(argv[1]); + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + { + PlasticityContactProblem<3> problem(prm); + problem.run(); + } + + return 0; } -- 2.39.5