From: frohne Date: Mon, 20 Aug 2012 17:53:21 +0000 (+0000) Subject: fix compiler warnings and renaming some classes and functions X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4c36a6fca4c2d7c72546bb74cd020cd3c1db2b35;p=dealii-svn.git fix compiler warnings and renaming some classes and functions git-svn-id: https://svn.dealii.org/trunk@26036 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 455cfae46a..281ed94ac9 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -63,1389 +63,1386 @@ #include #include - // This is new, however: in the previous - // example we got some unwanted output from - // the linear solvers. If we want to suppress - // it, we have to include this file and add a - // single line somewhere to the program (see - // the main() function below for that): + #include - // The final step, as in previous - // programs, is to import all the - // deal.II class and function names - // into the global namespace: -using namespace dealii; +namespace Step42 +{ + using namespace dealii; - // @sect3{The Step4 class template} + // @sect3{The PlasticityContactProblem class template} -template class ConstitutiveLaw; + template class ConstitutiveLaw; -template -class Step4 -{ -public: - Step4 (int _n_refinements_global, int _n_refinements_local); - void run (); - -private: - void make_grid (); - void setup_system(); - void assemble_mass_matrix (); - void assemble_nl_system (TrilinosWrappers::MPI::Vector &u); - void residual_nl_system (TrilinosWrappers::MPI::Vector &u, - Vector &sigma_eff_vector); - void projection_active_set (); - void dirichlet_constraints (); - void solve (); - void solve_newton (); - void output_results (const std::string& title) const; - void move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const; - void output_results (TrilinosWrappers::MPI::Vector vector, const std::string& title) const; - void output_results (Vector vector, const std::string& title) const; - - MPI_Comm mpi_communicator; - - parallel::distributed::Triangulation triangulation; - - FESystem fe; - DoFHandler dof_handler; - - IndexSet locally_owned_dofs; - IndexSet locally_relevant_dofs; - - int n_refinements_global; - int n_refinements_local; - unsigned int number_iterations; - std::vector run_time; - - ConstraintMatrix constraints; - ConstraintMatrix constraints_hanging_nodes; - ConstraintMatrix constraints_dirichlet_hanging_nodes; - - TrilinosWrappers::SparseMatrix system_matrix_newton; - TrilinosWrappers::SparseMatrix mass_matrix; - - TrilinosWrappers::MPI::Vector solution; - TrilinosWrappers::MPI::Vector old_solution; - TrilinosWrappers::MPI::Vector system_rhs_newton; - TrilinosWrappers::MPI::Vector resid_vector; - TrilinosWrappers::MPI::Vector diag_mass_matrix_vector; - IndexSet active_set; - - ConditionalOStream pcout; - - TrilinosWrappers::PreconditionAMG::AdditionalData additional_data; - TrilinosWrappers::PreconditionAMG preconditioner_u; - TrilinosWrappers::PreconditionAMG preconditioner_t; - - std::auto_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 - - std_cxx1x::shared_ptr Mp_preconditioner; -}; - -template -class ConstitutiveLaw -{ -public: - ConstitutiveLaw (double _E, double _nu, double _sigma_0, double _gamma, MPI_Comm _mpi_communicator, ConditionalOStream _pcout); - // ConstitutiveLaw (double mu, double kappa); - void plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, - SymmetricTensor<2,dim> &strain_tensor, - unsigned int &elast_points, - unsigned int &plast_points, - double &sigma_eff, - double &yield); - void linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, - SymmetricTensor<4,dim> &stress_strain_tensor, - 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; - -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 -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)); - pcout<< "-----> mu = " << mu << ", kappa = " << kappa <(), unit_symmetric_tensor()); - stress_strain_tensor_mu = 2*mu*(identity_tensor() - outer_product(unit_symmetric_tensor(), unit_symmetric_tensor())/3.0); -} + template + class PlasticityContactProblem + { + public: + PlasticityContactProblem (int _n_refinements_global, int _n_refinements_local); + void run (); -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; + private: + void make_grid (); + void setup_system(); + void assemble_mass_matrix (); + void assemble_nl_system (TrilinosWrappers::MPI::Vector &u); + void residual_nl_system (TrilinosWrappers::MPI::Vector &u, + Vector &sigma_eff_vector); + void update_solution_and_constraints (); + void dirichlet_constraints (); + void solve (); + void solve_newton (); + void output_results (const std::string& title) const; + void move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const; + void output_results (TrilinosWrappers::MPI::Vector vector, const std::string& title) const; + void output_results (Vector vector, const std::string& title) const; - tmp = fe_values[displacement].symmetric_gradient (shape_func,q_point); + int n_refinements_global; + int n_refinements_local; - return tmp; -} + MPI_Comm mpi_communicator; -template -void ConstitutiveLaw::plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, - SymmetricTensor<2,dim> &strain_tensor, - unsigned int &elast_points, - unsigned int &plast_points, - double &sigma_eff, - double &yield) -{ - if (dim == 3) - { - SymmetricTensor<2,dim> stress_tensor; - stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; - double tmp = E/((1+nu)*(1-2*nu)); + parallel::distributed::Triangulation triangulation; - SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); + FESystem fe; + DoFHandler dof_handler; - double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; - yield = 0; - stress_strain_tensor = stress_strain_tensor_mu; - double beta = 1.0; - if (deviator_stress_tensor_norm >= sigma_0) - { - beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; - stress_strain_tensor *= beta; - yield = 1; - plast_points += 1; - } - else - elast_points += 1; + unsigned int number_iterations; + std::vector run_time; -// std::cout<< beta < -void ConstitutiveLaw::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, - SymmetricTensor<4,dim> &stress_strain_tensor, - 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; - double tmp = E/((1+nu)*(1-2*nu)); + TrilinosWrappers::MPI::Vector solution; + TrilinosWrappers::MPI::Vector old_solution; + TrilinosWrappers::MPI::Vector system_rhs_newton; + TrilinosWrappers::MPI::Vector resid_vector; + TrilinosWrappers::MPI::Vector diag_mass_matrix_vector; + IndexSet active_set; - stress_strain_tensor = stress_strain_tensor_mu; - stress_strain_tensor_linearized = stress_strain_tensor_mu; + ConditionalOStream pcout; - SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); + TrilinosWrappers::PreconditionAMG::AdditionalData additional_data; + TrilinosWrappers::PreconditionAMG preconditioner_u; + TrilinosWrappers::PreconditionAMG preconditioner_t; - double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); + std::auto_ptr > plast_lin_hard; - double beta = 1.0; - if (deviator_stress_tensor_norm >= sigma_0) - { - beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; - stress_strain_tensor *= beta; - stress_strain_tensor_linearized *= beta; - deviator_stress_tensor /= deviator_stress_tensor_norm; - stress_strain_tensor_linearized -= beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); - } + double sigma_0; // Yield stress + double gamma; // Parameter for the linear isotropic hardening + double e_modul; // E-Modul + double nu; // Poisson ratio - stress_strain_tensor += stress_strain_tensor_kappa; - stress_strain_tensor_linearized += stress_strain_tensor_kappa; - } -} + std_cxx1x::shared_ptr Mp_preconditioner; + }; -namespace EquationData -{ template - class RightHandSide : public Function + class ConstitutiveLaw { 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; + ConstitutiveLaw (double _E, double _nu, double _sigma_0, double _gamma, MPI_Comm _mpi_communicator, ConditionalOStream _pcout); + // ConstitutiveLaw (double mu, double kappa); + void plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, + SymmetricTensor<2,dim> &strain_tensor, + unsigned int &elast_points, + unsigned int &plast_points, + double &sigma_eff, + double &yield); + void linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, + SymmetricTensor<4,dim> &stress_strain_tensor, + 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; + + 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 - double RightHandSide::value (const Point &p, - const unsigned int component) const + 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) { - double return_value = 0.0; - - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - // if ((p(0)-0.5)*(p(0)-0.5)+(p(1)-0.5)*(p(1)-0.5) < 0.2) - // return_value = -5000; - // else - return_value = 0.0; - // for (unsigned int i=0; i mu = " << mu << ", kappa = " << kappa <(), unit_symmetric_tensor()); + stress_strain_tensor_mu = 2*mu*(identity_tensor() - outer_product(unit_symmetric_tensor(), unit_symmetric_tensor())/3.0); } template - void RightHandSide::vector_value (const Point &p, - Vector &values) const + inline + SymmetricTensor<2,dim> ConstitutiveLaw::get_strain (const FEValues &fe_values, + const unsigned int shape_func, + const unsigned int q_point) const { - for (unsigned int c=0; cn_components; ++c) - values(c) = RightHandSide::value (p, c); - } + const FEValuesExtractors::Vector displacement (0); + SymmetricTensor<2,dim> tmp; + + tmp = fe_values[displacement].symmetric_gradient (shape_func,q_point); + return tmp; + } template - class BoundaryValues : public Function + void ConstitutiveLaw::plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, + SymmetricTensor<2,dim> &strain_tensor, + unsigned int &elast_points, + unsigned int &plast_points, + double &sigma_eff, + double &yield) { - public: - BoundaryValues () : Function(dim) {}; + if (dim == 3) + { + SymmetricTensor<2,dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; + double tmp = E/((1+nu)*(1-2*nu)); - virtual double value (const Point &p, - const unsigned int component = 0) const; + SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); - virtual void vector_value (const Point &p, - Vector &values) const; - }; + double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); - template - double BoundaryValues::value (const Point &p, - const unsigned int component) const - { - double return_value = 0; + yield = 0; + stress_strain_tensor = stress_strain_tensor_mu; + double beta = 1.0; + if (deviator_stress_tensor_norm >= sigma_0) + { + beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; + stress_strain_tensor *= beta; + yield = 1; + plast_points += 1; + } + else + elast_points += 1; - if (component == 0) - return_value = 0.0; - if (component == 1) - return_value = 0.0; - if (component == 2) - return_value = 0.0; + // std::cout<< beta < - void BoundaryValues::vector_value (const Point &p, - Vector &values) const + void ConstitutiveLaw::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, + SymmetricTensor<4,dim> &stress_strain_tensor, + SymmetricTensor<2,dim> &strain_tensor) { - for (unsigned int c=0; cn_components; ++c) - values(c) = BoundaryValues::value (p, c); - } + if (dim == 3) + { + SymmetricTensor<2,dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; + double tmp = E/((1+nu)*(1-2*nu)); + stress_strain_tensor = stress_strain_tensor_mu; + stress_strain_tensor_linearized = stress_strain_tensor_mu; - template - class Obstacle : public Function - { - public: - Obstacle () : Function(dim) {}; + SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &values) const; - }; + double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); - template - double Obstacle::value (const Point &p, - const unsigned int component) const - { - double R = 0.03; - double return_value = 0.0; - if (component == 0) - return_value = p(0); - if (component == 1) - return_value = p(1); - if (component == 2) + double beta = 1.0; + if (deviator_stress_tensor_norm >= sigma_0) { - // double hz = 0.98; - // double position_x = 0.5; - // double alpha = 12.0; - // double s_x = 0.5039649116; - // double s_y = hz + 0.00026316298; - // if (p(0) > position_x - R && p(0) < s_x) - // { - // return_value = -sqrt(R*R - (p(0)-position_x)*(p(0)-position_x)) + hz + R; - // } - // else if (p(0) >= s_x) - // { - // return_value = 12.0/90.0*p(0) + (s_y - alpha/90.0*s_x); - // } - // else - // return_value = 1e+10; - - // Hindernis Dortmund - double x1 = p(0); - double x2 = p(1); - if (((x2-0.5)*(x2-0.5)+(x1-0.5)*(x1-0.5)<=0.3*0.3)&&((x2-0.5)*(x2-0.5)+(x1-1.0)*(x1-1.0)>=0.4*0.4)&&((x2-0.5)*(x2-0.5)+x1*x1>=0.4*0.4)) - return_value = 0.999; - else - return_value = 1e+10; - - // Hindernis Werkzeug TKSE - // double shift_walze_x = 0.0; - // double shift_walze_y = 0.0; - // return_value = 0.032 + data->dicke - input_copy->mikro_height (p(0) + shift_walze_x, p(1) + shift_walze_y, p(2)); - - // Ball with radius R - // double R = 0.5; - // if (std::pow ((p(0)-1.0/2.0), 2) + std::pow ((p(1)-1.0/2.0), 2) < R*R) - // return_value = 1.0 + R - 0.001 - sqrt (R*R - std::pow ((p(0)-1.0/2.0), 2) - // - std::pow ((p(1)-1.0/2.0), 2)); - // else - // return_value = 1e+5; + beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; + stress_strain_tensor *= beta; + stress_strain_tensor_linearized *= beta; + deviator_stress_tensor /= deviator_stress_tensor_norm; + stress_strain_tensor_linearized -= beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); } - return return_value; - // return 1e+10;//0.98; + stress_strain_tensor += stress_strain_tensor_kappa; + stress_strain_tensor_linearized += stress_strain_tensor_kappa; + } } - template - void Obstacle::vector_value (const Point &p, - Vector &values) const + namespace EquationData { - for (unsigned int c=0; cn_components; ++c) - values(c) = Obstacle::value (p, c); - } -} + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim) {} + virtual double value (const Point &p, + const unsigned int component = 0) const; - // @sect3{Implementation of the Step4 class} - - // Next for the implementation of the class - // template that makes use of the functions - // above. As before, we will write everything - -template -Step4::Step4 (int _n_refinements_global, int _n_refinements_local) - : - n_refinements_global (_n_refinements_global), - n_refinements_local (_n_refinements_local), - mpi_communicator (MPI_COMM_WORLD), - triangulation (mpi_communicator), - fe (FE_Q(1), dim), - dof_handler (triangulation), - pcout (std::cout, - (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), - sigma_0 (400), - gamma (1.e-2), - e_modul (2.0e5), - nu (0.3) -{ - // double _E, double _nu, double _sigma_0, double _gamma - plast_lin_hard.reset (new ConstitutiveLaw (e_modul, nu, sigma_0, gamma, mpi_communicator, pcout)); -} + virtual void vector_value (const Point &p, + Vector &values) const; + }; -template -void Step4::make_grid () -{ - std::vector repet(3); - repet[0] = 1;//20; - repet[1] = 1; - repet[2] = 1; - - Point p1 (0,0,0); - Point p2 (1.0, 1.0, 1.0); - GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2); - - Triangulation<3>::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - - /* boundary_indicators: - _______ - / 9 /| - /______ / | - 8| | 8| - | 8 | / - |_______|/ - 6 - */ - - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - { - if (cell->face (face)->center ()[2] == p2(2)) - cell->face (face)->set_boundary_indicator (9); - 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); - } + 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) + // if ((p(0)-0.5)*(p(0)-0.5)+(p(1)-0.5)*(p(1)-0.5) < 0.2) + // return_value = -5000; + // else + return_value = 0.0; + // for (unsigned int i=0; i + void RightHandSide::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = RightHandSide::value (p, c); + } - triangulation.refine_global (n_refinements_global); - // Lokale Verfeinerung des Gitters - for (int step=0; step + class BoundaryValues : public Function { - cell = triangulation.begin_active(); // Iterator ueber alle Zellen - - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - { -// if (cell->face (face)->at_boundary() -// && cell->face (face)->boundary_indicator () == 9) -// { -// cell->set_refine_flag (); -// break; -// } -// else if (cell->level () == n_refinements + n_refinements_local - 1) -// { -// cell->set_refine_flag (); -// break; -// } - -// if (cell->face (face)->at_boundary() -// && cell->face (face)->boundary_indicator () == 9) -// { -// if (cell->face (face)->vertex (0)(0) <= 0.7 && -// cell->face (face)->vertex (1)(0) >= 0.3 && -// cell->face (face)->vertex (0)(1) <= 0.875 && -// cell->face (face)->vertex (2)(1) >= 0.125) -// { -// cell->set_refine_flag (); -// break; -// } -// } - - if (step == 0 && - cell->center ()(2) < n_refinements_local*9.0/64.0) - { - cell->set_refine_flag (); - break; - } - }; - triangulation.execute_coarsening_and_refinement (); + 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 -void Step4::setup_system () -{ - // setup dofs - { - dof_handler.distribute_dofs (fe); + template + double BoundaryValues::value (const Point &p, + const unsigned int component) const + { + double return_value = 0; - locally_owned_dofs = dof_handler.locally_owned_dofs (); - locally_relevant_dofs.clear(); - DoFTools::extract_locally_relevant_dofs (dof_handler, - locally_relevant_dofs); + 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; cn_components; ++c) + values(c) = BoundaryValues::value (p, c); + } + + + template + class Obstacle : public Function + { + public: + Obstacle () : 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 Obstacle::value (const Point &p, + const unsigned int component) const + { + double R = 0.03; + double return_value = 0.0; + if (component == 0) + return_value = p(0); + if (component == 1) + return_value = p(1); + if (component == 2) + { + // double hz = 0.98; + // double position_x = 0.5; + // double alpha = 12.0; + // double s_x = 0.5039649116; + // double s_y = hz + 0.00026316298; + // if (p(0) > position_x - R && p(0) < s_x) + // { + // return_value = -sqrt(R*R - (p(0)-position_x)*(p(0)-position_x)) + hz + R; + // } + // else if (p(0) >= s_x) + // { + // return_value = 12.0/90.0*p(0) + (s_y - alpha/90.0*s_x); + // } + // else + // return_value = 1e+10; + + // Hindernis Dortmund + double x1 = p(0); + double x2 = p(1); + if (((x2-0.5)*(x2-0.5)+(x1-0.5)*(x1-0.5)<=0.3*0.3)&&((x2-0.5)*(x2-0.5)+(x1-1.0)*(x1-1.0)>=0.4*0.4)&&((x2-0.5)*(x2-0.5)+x1*x1>=0.4*0.4)) + return_value = 0.999; + else + return_value = 1e+10; + + // Hindernis Werkzeug TKSE + // double shift_walze_x = 0.0; + // double shift_walze_y = 0.0; + // return_value = 0.032 + data->dicke - input_copy->mikro_height (p(0) + shift_walze_x, p(1) + shift_walze_y, p(2)); + + // Ball with radius R + // double R = 0.5; + // if (std::pow ((p(0)-1.0/2.0), 2) + std::pow ((p(1)-1.0/2.0), 2) < R*R) + // return_value = 1.0 + R - 0.001 - sqrt (R*R - std::pow ((p(0)-1.0/2.0), 2) + // - std::pow ((p(1)-1.0/2.0), 2)); + // else + // return_value = 1e+5; + } + return return_value; + + // return 1e+10;//0.98; + } + + template + void Obstacle::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = Obstacle::value (p, c); + } } - // setup hanging nodes and dirichlet constraints + + // @sect3{Implementation of the PlasticityContactProblem class} + + // Next for the implementation of the class + // template that makes use of the functions + // above. As before, we will write everything + + template + PlasticityContactProblem::PlasticityContactProblem (int _n_refinements_global, int _n_refinements_local) + : + n_refinements_global (_n_refinements_global), + n_refinements_local (_n_refinements_local), + mpi_communicator (MPI_COMM_WORLD), + triangulation (mpi_communicator), + fe (FE_Q(1), dim), + dof_handler (triangulation), + pcout (std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), + sigma_0 (400), + gamma (1.e-2), + e_modul (2.0e5), + nu (0.3) { - // constraints_hanging_nodes.clear (); - 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_active_cells() - << std::endl - << "Total number of cells: " - << triangulation.n_cells() - << std::endl - << "Number of degrees of freedom: " - << dof_handler.n_dofs () - << std::endl; - - dirichlet_constraints (); + // double _E, double _nu, double _sigma_0, double _gamma + plast_lin_hard.reset (new ConstitutiveLaw (e_modul, nu, sigma_0, gamma, mpi_communicator, pcout)); } - // Initialzation for matrices and vectors + template + void PlasticityContactProblem::make_grid () { - solution.reinit (locally_relevant_dofs, mpi_communicator); - system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator); - old_solution.reinit (system_rhs_newton); - resid_vector.reinit (system_rhs_newton); - diag_mass_matrix_vector.reinit (system_rhs_newton); - active_set.set_size (locally_relevant_dofs.size ()); + std::vector repet(3); + repet[0] = 1;//20; + repet[1] = 1; + repet[2] = 1; + + Point p1 (0,0,0); + Point p2 (1.0, 1.0, 1.0); + GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2); + + Triangulation<3>::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + + /* boundary_indicators: + _______ + / 9 /| + /______ / | + 8| | 8| + | 8 | / + |_______|/ + 6 + */ + + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face (face)->center ()[2] == p2(2)) + cell->face (face)->set_boundary_indicator (9); + 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_refinements_global); + + // Lokale Verfeinerung des Gitters + for (int step=0; step::faces_per_cell; ++face) + { + // if (cell->face (face)->at_boundary() + // && cell->face (face)->boundary_indicator () == 9) + // { + // cell->set_refine_flag (); + // break; + // } + // else if (cell->level () == n_refinements + n_refinements_local - 1) + // { + // cell->set_refine_flag (); + // break; + // } + + // if (cell->face (face)->at_boundary() + // && cell->face (face)->boundary_indicator () == 9) + // { + // if (cell->face (face)->vertex (0)(0) <= 0.7 && + // cell->face (face)->vertex (1)(0) >= 0.3 && + // cell->face (face)->vertex (0)(1) <= 0.875 && + // cell->face (face)->vertex (2)(1) >= 0.125) + // { + // cell->set_refine_flag (); + // break; + // } + // } + + if (step == 0 && + cell->center ()(2) < n_refinements_local*9.0/64.0) + { + cell->set_refine_flag (); + break; + } + }; + triangulation.execute_coarsening_and_refinement (); + }; } - // setup sparsity pattern + template + void PlasticityContactProblem::setup_system () { - TrilinosWrappers::SparsityPattern sp (locally_owned_dofs, - mpi_communicator); + // setup dofs + { + dof_handler.distribute_dofs (fe); - DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false, - Utilities::MPI::this_mpi_process(mpi_communicator)); + locally_owned_dofs = dof_handler.locally_owned_dofs (); + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + } - sp.compress(); + // setup hanging nodes and dirichlet constraints + { + // constraints_hanging_nodes.clear (); + 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_active_cells() + << std::endl + << "Total number of cells: " + << triangulation.n_cells() + << std::endl + << "Number of degrees of freedom: " + << dof_handler.n_dofs () + << std::endl; + + dirichlet_constraints (); + } - system_matrix_newton.reinit (sp); + // Initialzation for matrices and vectors + { + solution.reinit (locally_relevant_dofs, mpi_communicator); + system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator); + old_solution.reinit (system_rhs_newton); + resid_vector.reinit (system_rhs_newton); + diag_mass_matrix_vector.reinit (system_rhs_newton); + active_set.set_size (locally_relevant_dofs.size ()); + } - mass_matrix.reinit (sp); - } + // setup sparsity pattern + { + TrilinosWrappers::SparsityPattern sp (locally_owned_dofs, + mpi_communicator); - assemble_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 -void Step4::assemble_mass_matrix () -{ - QTrapez face_quadrature_formula; + sp.compress(); + + system_matrix_newton.reinit (sp); - FEFaceValues fe_values_face (fe, face_quadrature_formula, - update_values | update_quadrature_points | update_JxW_values); + mass_matrix.reinit (sp); + } - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int dofs_per_face = fe.dofs_per_face; - const unsigned int n_face_q_points = face_quadrature_formula.size(); + assemble_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 cell_matrix (dofs_per_cell, dofs_per_cell); + diag_mass_matrix_vector.compress (); + } - std::vector local_dof_indices (dofs_per_cell); + template + void PlasticityContactProblem::assemble_mass_matrix () + { + QTrapez face_quadrature_formula; - const FEValuesExtractors::Vector displacement (0); + FEFaceValues fe_values_face (fe, face_quadrature_formula, + update_values | update_quadrature_points | update_JxW_values); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int dofs_per_face = fe.dofs_per_face; + const unsigned int n_face_q_points = face_quadrature_formula.size(); - for (; cell!=endc; ++cell) - if (cell->is_locally_owned()) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face (face)->at_boundary() - && cell->face (face)->boundary_indicator () == 9) - { - fe_values_face.reinit (cell, face); - cell_matrix = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix, - local_dof_indices, - mass_matrix); - } - - mass_matrix.compress (); -} + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); -template -void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) -{ - QGauss quadrature_formula(2); - QGauss face_quadrature_formula(2); + std::vector local_dof_indices (dofs_per_cell); - FEValues fe_values (fe, quadrature_formula, - UpdateFlags(update_values | - update_gradients | - update_q_points | - update_JxW_values)); + const FEValuesExtractors::Vector displacement (0); - FEFaceValues fe_values_face (fe, face_quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size (); - const unsigned int n_face_q_points = face_quadrature_formula.size(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face (face)->at_boundary() + && cell->face (face)->boundary_indicator () == 9) + { + fe_values_face.reinit (cell, face); + cell_matrix = 0; - const EquationData::RightHandSide right_hand_side; - std::vector > right_hand_side_values (n_q_points, - Vector(dim)); - std::vector > right_hand_side_values_face (n_face_q_points, - Vector(dim)); + for (unsigned int q_point=0; q_point cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); - std::vector local_dof_indices (dofs_per_cell); + constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix, + local_dof_indices, + mass_matrix); + } - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + mass_matrix.compress (); + } - const FEValuesExtractors::Vector displacement (0); + template + void PlasticityContactProblem::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) + { + QGauss quadrature_formula(2); + QGauss face_quadrature_formula(2); + + 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); + + TrilinosWrappers::MPI::Vector test_rhs(solution); + 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 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]); + + // if (q_point == 0) + // std::cout<< stress_strain_tensor_linearized <get_strain(fe_values, i, q_point); + + for (unsigned int j=0; jget_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::faces_per_cell; ++face) + { + if (cell->face (face)->at_boundary() + && cell->face (face)->boundary_indicator () == 9) + { + 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 rhs_values; + rhs_values = 0; + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix_newton, system_rhs_newton, true); + }; - TrilinosWrappers::MPI::Vector test_rhs(solution); - 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 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]); - - // if (q_point == 0) - // std::cout<< stress_strain_tensor_linearized <get_strain(fe_values, i, q_point); - - for (unsigned int j=0; jget_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::faces_per_cell; ++face) - { - if (cell->face (face)->at_boundary() - && cell->face (face)->boundary_indicator () == 9) - { - 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 rhs_values; - rhs_values = 0; - for (unsigned int i=0; iget_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 (); + system_rhs_newton.compress (Add); + } - system_matrix_newton.compress (); - system_rhs_newton.compress (Add); -} + template + void PlasticityContactProblem::residual_nl_system (TrilinosWrappers::MPI::Vector &u, + Vector &sigma_eff_vector) + { + QGauss quadrature_formula(2); + QGauss face_quadrature_formula(2); + + 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); + Vector cell_sigma_eff (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 sigma_eff = 0; + double yield = 0; + unsigned int cell_number = 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 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, sigma_eff, yield); + + // sigma_eff_vector (cell_number) += sigma_eff; + sigma_eff_vector (cell_number) += yield; + + /* if (q_point == 0) + std::cout<< stress_strain_tensor <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::faces_per_cell; ++face) + { + if (cell->face (face)->at_boundary() + && cell->face (face)->boundary_indicator () == 9) + { + 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 rhs_values; + rhs_values = 0; + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs, + local_dof_indices, + system_rhs_newton); + + sigma_eff_vector(cell_number) /= n_q_points; + cell_number += 1; + }; -template -void Step4::residual_nl_system (TrilinosWrappers::MPI::Vector &u, - Vector &sigma_eff_vector) -{ - QGauss quadrature_formula(2); - QGauss face_quadrature_formula(2); - - 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); - Vector cell_sigma_eff (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 sigma_eff = 0; - double yield = 0; - unsigned int cell_number = 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 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, sigma_eff, yield); - - // sigma_eff_vector (cell_number) += sigma_eff; - sigma_eff_vector (cell_number) += yield; - - /* if (q_point == 0) - std::cout<< stress_strain_tensor <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::faces_per_cell; ++face) - { - if (cell->face (face)->at_boundary() - && cell->face (face)->boundary_indicator () == 9) - { - 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 rhs_values; - rhs_values = 0; - for (unsigned int i=0; iget_dof_indices (local_dof_indices); - constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs, - local_dof_indices, - system_rhs_newton); - - sigma_eff_vector(cell_number) /= n_q_points; - cell_number += 1; - }; + system_rhs_newton.compress (); - system_rhs_newton.compress (); + unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator); + unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator); + pcout<< "Elast-Points = " << sum_elast_points < obstacle; + std::vector vertex_touched (dof_handler.n_dofs (), false); -template -void Step4::dirichlet_constraints () -{ - /* boundary_indicators: - _______ - / 9 /| - /______ / | - 8| | 8| - | 8 | / - |_______|/ - 6 - */ - - constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs); - constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes); - - std::vector component_mask (dim, true); - component_mask[0] = true; - component_mask[1] = true; - component_mask[2] = true; - VectorTools::interpolate_boundary_values (dof_handler, - 6, - EquationData::BoundaryValues(), - constraints_dirichlet_hanging_nodes, - component_mask); - - component_mask[0] = true; - component_mask[1] = true; - component_mask[2] = false; - VectorTools::interpolate_boundary_values (dof_handler, - 8, - EquationData::BoundaryValues(), - constraints_dirichlet_hanging_nodes, - component_mask); - constraints_dirichlet_hanging_nodes.close (); -} + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); -template -void Step4::solve () -{ - pcout << "Solving ..." << std::endl; - Timer t; + 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; + + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face (face)->at_boundary() + && cell->face (face)->boundary_indicator () == 9) + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + unsigned int index_z = cell->face (face)->vertex_dof_index (v,2); + + if (vertex_touched[cell->face (face)->vertex_index(v)] == false) + vertex_touched[cell->face (face)->vertex_index(v)] = true; + else + continue; + + // the local row where + Point point (cell->face (face)->vertex (v)[0],/* + solution (index_x),*/ + cell->face (face)->vertex (v)[1], + cell->face (face)->vertex (v)[2]); + + double obstacle_value = obstacle.value (point, 2); + double solution_index_z = solution (index_z); + double gap = obstacle_value - point (2); + + if (lambda (index_z) + + c * + diag_mass_matrix_vector_relevant (index_z) * + (solution_index_z - gap) + > 0) + { + constraints.add_line (index_z); + constraints.set_inhomogeneity (index_z, gap); + + distributed_solution (index_z) = gap; + + if (locally_relevant_dofs.is_element (index_z)) + active_set.add_index (index_z); + + if (locally_owned_dofs.is_element (index_z)) + active_set_locally_owned.add_index (index_z); + + // std::cout<< index_z << ", " + // << "Error: " << lambda (index_z) + + // diag_mass_matrix_vector_relevant (index_z)*c*(solution_index_z - gap) + // << ", " << lambda (index_z) + // << ", " << diag_mass_matrix_vector_relevant (index_z) + // << ", " << obstacle_value + // << ", " << solution_index_z + // < + void PlasticityContactProblem::dirichlet_constraints () + { + /* boundary_indicators: + _______ + / 9 /| + /______ / | + 8| | 8| + | 8 | / + |_______|/ + 6 + */ + + constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs); + constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes); + + std::vector component_mask (dim, true); + component_mask[0] = true; + component_mask[1] = true; + component_mask[2] = true; + VectorTools::interpolate_boundary_values (dof_handler, + 6, + EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, + component_mask); + + component_mask[0] = true; + component_mask[1] = true; + component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 8, + EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, + component_mask); + constraints_dirichlet_hanging_nodes.close (); + } - // Solving iterative + template + void PlasticityContactProblem::solve () + { + pcout << "Solving ..." << std::endl; + Timer t; - MPI_Barrier (mpi_communicator); - t.restart(); + TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); + distributed_solution = solution; - preconditioner_u.initialize (system_matrix_newton, additional_data); - - MPI_Barrier (mpi_communicator); - t.stop(); - if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) - run_time[6] += t.wall_time(); - - MPI_Barrier (mpi_communicator); - t.restart(); - -// ReductionControl reduction_control (10000, 1e-15, 1e-4); -// SolverCG -// solver (reduction_control, mpi_communicator); -// solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); - - PrimitiveVectorMemory mem; - TrilinosWrappers::MPI::Vector tmp (system_rhs_newton); - const double solver_tolerance = 1e-4 * - system_matrix_newton.residual (tmp, distributed_solution, system_rhs_newton); - SolverControl solver_control (system_matrix_newton.m(), solver_tolerance); - SolverFGMRES - solver(solver_control, mem, - SolverFGMRES:: - AdditionalData(30, true)); - solver.solve(system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); - - pcout << "Initial error: " << solver_control.initial_value() < -void Step4::solve_newton () -{ - double resid=0; - double resid_old=100000; - TrilinosWrappers::MPI::Vector res (system_rhs_newton); - TrilinosWrappers::MPI::Vector tmp_vector (system_rhs_newton); - Timer t; - - std::vector > constant_modes; - std::vector components (dim,true); - components[dim] = false; - DoFTools::extract_constant_modes (dof_handler, components, - 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); - Vector sigma_eff_vector; - sigma_eff_vector.reinit (triangulation.n_active_cells()); - unsigned int j = 0; - unsigned int number_assemble_system = 0; - for (; j<=100;j++) - { - pcout<< " " <(i)); - old_solution = tmp_vector; - old_solution.sadd(1-a,a, distributed_solution); - - MPI_Barrier (mpi_communicator); - t.restart(); - system_rhs_newton = 0; - sigma_eff_vector = 0; - solution = old_solution; - residual_nl_system (solution, sigma_eff_vector); - 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::type_dof_data, - data_component_interpretation); - data_out.add_data_vector (lambda, std::vector(dim, "Residual"), - DataOut::type_dof_data, - data_component_interpretation); - data_out.add_data_vector (active_set, std::vector(dim, "ActiveSet"), - DataOut::type_dof_data, - data_component_interpretation); + PrimitiveVectorMemory mem; + TrilinosWrappers::MPI::Vector tmp (system_rhs_newton); + const double solver_tolerance = 1e-4 * + system_matrix_newton.residual (tmp, distributed_solution, system_rhs_newton); + SolverControl solver_control (system_matrix_newton.m(), solver_tolerance); + SolverFGMRES + solver(solver_control, mem, + SolverFGMRES:: + AdditionalData(30, true)); + solver.solve(system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); - Vector subdomain (triangulation.n_active_cells()); - for (unsigned int i=0; i filenames; - for (unsigned int i=0; - i + void PlasticityContactProblem::solve_newton () + { + double resid=0; + double resid_old=100000; + TrilinosWrappers::MPI::Vector res (system_rhs_newton); + TrilinosWrappers::MPI::Vector tmp_vector (system_rhs_newton); + Timer t; + + std::vector > constant_modes; + std::vector components (dim,true); + components[dim] = false; + DoFTools::extract_constant_modes (dof_handler, components, + 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); + Vector sigma_eff_vector; + sigma_eff_vector.reinit (triangulation.n_active_cells()); + unsigned int j = 0; + unsigned int number_assemble_system = 0; + for (; j<=100;j++) + { + pcout<< " " <(i)); + old_solution = tmp_vector; + old_solution.sadd(1-a,a, distributed_solution); + + MPI_Barrier (mpi_communicator); + t.restart(); + system_rhs_newton = 0; + sigma_eff_vector = 0; + solution = old_solution; + residual_nl_system (solution, sigma_eff_vector); + 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 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::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; dvertex_dof_index(v,d)) != 0) - vertex_displacement[d] - = _complete_displacement(cell->vertex_dof_index(v,d)); - } - - cell->vertex(v) += vertex_displacement; - } - } -} + template + void PlasticityContactProblem::output_results (const std::string& title) const + { + move_mesh (solution); -template -void Step4::output_results (TrilinosWrappers::MPI::Vector vector, - const std::string& title) const -{ - DataOut data_out; + TrilinosWrappers::MPI::Vector lambda (solution); + lambda = resid_vector; - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (vector, "vector_to_plot"); + DataOut data_out; - data_out.build_patches (); + data_out.attach_dof_handler (dof_handler); - std::ofstream output_vtk (dim == 2 ? - (title + ".vtk").c_str () : - (title + ".vtk").c_str ()); - data_out.write_vtk (output_vtk); -} + const std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_out.add_data_vector (solution, std::vector(dim, "Displacement"), + DataOut::type_dof_data, + data_component_interpretation); + data_out.add_data_vector (lambda, std::vector(dim, "Residual"), + DataOut::type_dof_data, + data_component_interpretation); + data_out.add_data_vector (active_set, std::vector(dim, "ActiveSet"), + DataOut::type_dof_data, + data_component_interpretation); -template -void Step4::output_results (Vector vector, const std::string& title) const -{ - DataOut data_out; + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i -void Step4::run () -{ - pcout << "Solving problem in " << dim << " space dimensions." << std::endl; + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + std::vector filenames; + for (unsigned int i=0; + i + void PlasticityContactProblem::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const + { + pcout<< "Moving mesh." < 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::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; dvertex_dof_index(v,d)) != 0) + vertex_displacement[d] + = _complete_displacement(cell->vertex_dof_index(v,d)); + } + + cell->vertex(v) += vertex_displacement; + } + } + } - start = clock(); - make_grid(); - // mesh_surface (); + template + void PlasticityContactProblem::output_results (TrilinosWrappers::MPI::Vector vector, + const std::string& title) const + { + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (vector, "vector_to_plot"); + + data_out.build_patches (); + + std::ofstream output_vtk (dim == 2 ? + (title + ".vtk").c_str () : + (title + ".vtk").c_str ()); + data_out.write_vtk (output_vtk); + } + + template + void PlasticityContactProblem::output_results (Vector vector, const std::string& title) const + { + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (vector, "vector_to_plot"); + + data_out.build_patches (); + + std::ofstream output_vtk (dim == 2 ? + (title + ".vtk").c_str () : + (title + ".vtk").c_str ()); + data_out.write_vtk (output_vtk); + } + + template + void PlasticityContactProblem::run () + { + pcout << "Solving problem in " << dim << " space dimensions." << std::endl; - setup_system (); + run_time.resize (8); - end = clock(); - run_time[0] = (double)(end-start)/CLOCKS_PER_SEC; + clock_t start, end; + + start = clock(); + make_grid(); + // mesh_surface (); - solve_newton (); -} + setup_system (); + + end = clock(); + run_time[0] = (double)(end-start)/CLOCKS_PER_SEC; + solve_newton (); + } +} // @sect3{The main function} - // And this is the main function. It also - // looks mostly like in step-3, but if you - // look at the code below, note how we first - // create a variable of type - // Step4@<2@> (forcing - // the compiler to compile the class template - // with dim replaced by - // 2) and run a 2d simulation, - // and then we do the whole thing over in 3d. - // - // In practice, this is probably not what you - // would do very frequently (you probably - // either want to solve a 2d problem, or one - // in 3d, but not both at the same - // time). However, it demonstrates the - // mechanism by which we can simply change - // which dimension we want in a single place, - // and thereby force the compiler to - // recompile the dimension independent class - // templates for the dimension we - // request. The emphasis here lies on the - // fact that we only need to change a single - // place. This makes it rather trivial to - // debug the program in 2d where computations - // are fast, and then switch a single place - // to a 3 to run the much more computing - // intensive program in 3d for `real' - // computations. - // - // Each of the two blocks is enclosed in - // braces to make sure that the - // laplace_problem_2d variable - // goes out of scope (and releases the memory - // it holds) before we move on to allocate - // memory for the 3d case. Without the - // additional braces, the - // laplace_problem_2d variable - // would only be destroyed at the end of the - // function, i.e. after running the 3d - // problem, and would needlessly hog memory - // while the 3d run could actually use it. + // And this is the main function. It also + // looks mostly like in step-3, but if you + // look at the code below, note how we first + // create a variable of type + // PlasticityContactProblem@<2@> (forcing + // the compiler to compile the class template + // with dim replaced by + // 2) and run a 2d simulation, + // and then we do the whole thing over in 3d. + // + // In practice, this is probably not what you + // would do very frequently (you probably + // either want to solve a 2d problem, or one + // in 3d, but not both at the same + // time). However, it demonstrates the + // mechanism by which we can simply change + // which dimension we want in a single place, + // and thereby force the compiler to + // recompile the dimension independent class + // templates for the dimension we + // request. The emphasis here lies on the + // fact that we only need to change a single + // place. This makes it rather trivial to + // debug the program in 2d where computations + // are fast, and then switch a single place + // to a 3 to run the much more computing + // intensive program in 3d for `real' + // computations. + // + // Each of the two blocks is enclosed in + // braces to make sure that the + // laplace_problem_2d variable + // goes out of scope (and releases the memory + // it holds) before we move on to allocate + // memory for the 3d case. Without the + // additional braces, the + // laplace_problem_2d variable + // would only be destroyed at the end of the + // function, i.e. after running the 3d + // problem, and would needlessly hog memory + // while the 3d run could actually use it. // // Finally, the first line of the function is // used to suppress some output. Remember @@ -1480,6 +1477,9 @@ void Step4::run () // library. int main (int argc, char *argv[]) { + using namespace dealii; + using namespace Step42; + deallog.depth_console (0); clock_t start, end; @@ -1488,16 +1488,16 @@ int main (int argc, char *argv[]) Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv); { - int _n_refinements_global = 1; + int _n_refinements_global = 4; int _n_refinements_local = 1; - + if (argc == 3) { _n_refinements_global = atoi(argv[1]); _n_refinements_local = atoi(argv[2]); } - - Step4<3> laplace_problem_3d (_n_refinements_global, _n_refinements_local); + + PlasticityContactProblem<3> laplace_problem_3d (_n_refinements_global, _n_refinements_local); laplace_problem_3d.run (); }