From 242105b8fbdf5d6deac2f428b4484b8901264b17 Mon Sep 17 00:00:00 2001 From: frohne Date: Mon, 2 Apr 2012 13:00:37 +0000 Subject: [PATCH] Initial version of elasto-plastic contact problem git-svn-id: https://svn.dealii.org/trunk@25356 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/Makefile | 27 +- deal.II/examples/step-42/step-42.cc | 2424 +++++++++++++++------------ 2 files changed, 1392 insertions(+), 1059 deletions(-) diff --git a/deal.II/examples/step-42/Makefile b/deal.II/examples/step-42/Makefile index 8ffa568be9..4f520f02a3 100644 --- a/deal.II/examples/step-42/Makefile +++ b/deal.II/examples/step-42/Makefile @@ -30,7 +30,7 @@ D = ../../ # shall be deleted when calling `make clean'. Object and backup files, # executables and the like are removed anyway. Here, we give a list of # files in the various output formats that deal.II supports. -clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2 +clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *vtu *visit *ucd *.d2 @@ -45,25 +45,27 @@ clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2 # settings include $D/common/Make.global_options -################################################################ -# This example program will only work if UMFPACK is installed. If this +################################################################ +# This example program will only work if Trilinos is installed. If this # is not the case, then simply redefine the main targets to do nothing - - -ifneq ($(USE_CONTRIB_UMFPACK),yes) +ifneq ($(USE_CONTRIB_TRILINOS)$(USE_CONTRIB_P4EST),yesyes) default run clean: @echo - @echo "===========================================================" - @echo "= This program cannot be compiled without UMFPACK. Make =" - @echo "= sure you have UMFPACK installed and detected during =" - @echo "= configuration of deal.II =" + @echo "===========================================================" + @echo "= This program can only be run with the following =" + @echo "= configuration: =" + @echo "= - with Trilinos =" + @echo "= - with P4EST =" + @echo "= At least one of these conditions is not satisfied. =" @echo "===========================================================" @echo -endif +else # ################################################################ + + # Since the whole project consists of only one file, we need not # consider difficult dependencies. We only have to declare the # libraries which we want to link to the object file. deal.II has two @@ -108,6 +110,7 @@ run: $(target)$(EXEEXT) @./$(target)$(EXEEXT) + # As a last rule to the `make' program, we define what to do when # cleaning up a directory. This usually involves deleting object files # and other automatically created files such as the executable itself, @@ -130,6 +133,7 @@ clean: @$(CXX) $(CXXFLAGS.o) -c $< -o $@ + # The following statement tells make that the rules `run' and `clean' # are not expected to produce files of the same name as Makefile rules # usually do. @@ -160,3 +164,4 @@ Makefile.dep: $(target).cc Makefile \ include Makefile.dep +endif # USE_CONTRIB_TRILINOS diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 0c5fab4170..c4e9230ad1 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -1,1188 +1,1516 @@ -/* Author: Wolfgang Bangerth, Texas A&M University, 2008 */ + +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */ /* $Id$ */ /* */ -/* Copyright (C) 2008, 2009, 2010, 2011, 2012 by the deal.II authors */ +/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ -/* without copyright and license information. Please refer */ +/* without copyrightG and license information. Please refer */ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ - // @sect3{Include files} - // As usual, we start by including - // some well-known files: -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include + // The first few (many?) include + // files have already been used in + // the previous example, so we will + // not explain their meaning here + // again. #include +#include #include #include #include -#include -#include -#include -#include - -#include -#include #include -#include - +#include #include #include +#include #include -#include - +#include +#include #include #include -#include -#include - - -#include +#include +#include +#include +#include +#include +#include +#include +#include -#include +#include +#include +#include +#include +#include +#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include #include -#include - - -namespace Step42 -{ - using namespace dealii; +#include +#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; - template - struct InnerPreconditioner; + // @sect3{The Step4 class template} +template class ConstitutiveLaw; - template <> - struct InnerPreconditioner<2> - { - typedef SparseDirectUMFPACK type; - }; +template +class Step4 +{ +public: + Step4 (); + 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; + 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 +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 <> - struct InnerPreconditioner<3> - { - typedef SparseILU type; - }; + tmp = fe_values[displacement].symmetric_gradient (shape_func,q_point); + + return tmp; +} - template - void copy(const MATRIX &matrix, - FullMatrix &full_matrix) +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) +{ + // Plane strain + if (dim == 3) { - const unsigned int m = matrix.m(); - const unsigned int n = matrix.n(); - full_matrix.reinit(n,m); + 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)); + double stress_tensor_33 = 0.0;//tmp*(strain_tensor[0][0] + strain_tensor[1][1])*nu; + + SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); + + double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); + deviator_stress_tensor_norm = std::sqrt (deviator_stress_tensor_norm*deviator_stress_tensor_norm + + stress_tensor_33*stress_tensor_33); + + 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; + +// std::cout<< beta < unit (n); - Vector result (m); - for(unsigned int i=0; i - class StokesProblem +template +void ConstitutiveLaw::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, + SymmetricTensor<4,dim> &stress_strain_tensor, + SymmetricTensor<2,dim> &strain_tensor) +{ + // Plane strains + if (dim == 3) { - public: - StokesProblem (const unsigned int degree); - void run (); - - private: - void setup_dofs (); - void assemble_system (); - void assemble_multigrid (); - void solve (); - void solve_block (); - - void find_dofs_on_lower_level (std::vector > &lower_dofs, - std::vector > &boundary_dofs); - - void output_results (const unsigned int refinement_cycle) const; - void refine_mesh (); - - const unsigned int degree; - - Triangulation triangulation; - FESystem fe; - MGDoFHandler dof_handler; - - ConstraintMatrix constraints; - - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; - - BlockVector solution; - BlockVector system_rhs; - - MGLevelObject mg_constraints; - MGLevelObject mg_sparsity; - MGLevelObject > mg_matrices; - - MGLevelObject > mg_interface_matrices; - MGConstrainedDoFs mg_constrained_dofs; - std::vector > mg_dofs_per_component; - - std::vector::type> > mg_A_preconditioner; - std_cxx1x::shared_ptr::type> A_preconditioner; - }; - + 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)); + double stress_tensor_33 = 0.0;//tmp*(strain_tensor[0][0] + strain_tensor[1][1])*nu; + + stress_strain_tensor = stress_strain_tensor_mu; + stress_strain_tensor_linearized = stress_strain_tensor_mu; + + SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); + + double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); + deviator_stress_tensor_norm = std::sqrt (deviator_stress_tensor_norm*deviator_stress_tensor_norm + stress_tensor_33*stress_tensor_33); + 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); + } + + stress_strain_tensor += stress_strain_tensor_kappa; + stress_strain_tensor_linearized += stress_strain_tensor_kappa; + } +} +namespace EquationData +{ template - class BoundaryValues : public Function + class RightHandSide : public Function { - public: - BoundaryValues () : Function(dim+1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; + 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 - BoundaryValues::value (const Point &p, - const unsigned int component) const + double RightHandSide::value (const Point &p, + const unsigned int component) const { - Assert (component < this->n_components, - ExcIndexRange (component, 0, this->n_components)); - - if (component == 0 && p[0] == 0) - return (dim == 2 ? - p[1]*(p[1]-1.) : p[1]*(p[1]-1.) * p[2]*(p[2]-1.)); - return 0; + 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 - BoundaryValues::vector_value (const Point &p, - Vector &values) const + void RightHandSide::vector_value (const Point &p, + Vector &values) const { for (unsigned int c=0; cn_components; ++c) - values(c) = BoundaryValues::value (p, c); + values(c) = RightHandSide::value (p, c); } - - template - class RightHandSide : public Function + class BoundaryValues : public Function { - public: - RightHandSide () : Function(dim+1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; - + 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 - RightHandSide::value (const Point &/*p*/, - const unsigned int component) const + double BoundaryValues::value (const Point &p, + const unsigned int component) const { - return (component == 1 ? 1 : 0); - } + 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 - RightHandSide::vector_value (const Point &p, - Vector &values) const + void BoundaryValues::vector_value (const Point &p, + Vector &values) const { for (unsigned int c=0; cn_components; ++c) - values(c) = RightHandSide::value (p, c); + values(c) = BoundaryValues::value (p, c); } - - - - template - class InverseMatrix : public Subscriptor + + template + class Obstacle : public Function { - public: - InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner); + public: + Obstacle () : Function(dim) {}; - void vmult (Vector &dst, - const Vector &src) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - mutable std::string name; - private: - const SmartPointer matrix; - const SmartPointer preconditioner; + virtual void vector_value (const Point &p, + Vector &values) const; }; - - template - InverseMatrix::InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner) - : - matrix (&m), - preconditioner (&preconditioner) - {} - - - template - void InverseMatrix::vmult (Vector &dst, - const Vector &src) const + template + double Obstacle::value (const Point &p, + const unsigned int component) const { - SolverControl solver_control (src.size(), 1.0e-12*src.l2_norm()); - SolverCG<> cg (solver_control); - - dst = 0; - - try + 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) { - cg.solve (*matrix, dst, src, *preconditioner); + // 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.2; + 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; } - catch (...) - { - std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; - abort (); - } - -#ifdef STEP_42_TEST - if (name == "in schur") - std::cout << " " << solver_control.last_step() - << " inner CG steps inside the Schur complement "; - else if (name == "top left") - std::cout << " " << solver_control.last_step() - << " CG steps on the top left block "; - else if (name == "rhs") - std::cout << " " << solver_control.last_step() - << " CG steps for computing the r.h.s. "; - else - abort (); - - std::cout << solver_control.initial_value() << "->" << solver_control.last_value() - << std::endl; -#endif + return return_value; + + // return 1e+10;//0.98; } - - - template - class BlockSchurPreconditioner : public Subscriptor - { - public: - BlockSchurPreconditioner (const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner); - - void vmult (BlockVector &dst, - const BlockVector &src) const; - - private: - const SmartPointer > system_matrix; - const SmartPointer, - PreconditionerMp > > m_inverse; - const PreconditionerA &a_preconditioner; - - mutable Vector tmp; - - }; - - template - BlockSchurPreconditioner::BlockSchurPreconditioner( - const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner - ) - : - system_matrix (&S), - m_inverse (&Mpinv), - a_preconditioner (Apreconditioner), - tmp (S.block(1,1).m()) - {} - - // Now the interesting function, the multiplication of - // the preconditioner with a BlockVector. - template - void BlockSchurPreconditioner::vmult ( - BlockVector &dst, - const BlockVector &src) const - { - // Form u_new = A^{-1} u - a_preconditioner.vmult (dst.block(0), src.block(0)); - // Form tmp = - B u_new + p - // (SparseMatrix::residual - // does precisely this) - system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); - // Change sign in tmp - tmp *= -1; - // Multiply by approximate Schur complement - // (i.e. a pressure mass matrix) - m_inverse->vmult (dst.block(1), tmp); - } - - template - class SchurComplement : public Subscriptor - { - public: - SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix, Preconditioner> &A_inverse); - - void vmult (Vector &dst, - const Vector &src) const; - - unsigned int m() const - { - return system_matrix->block(1,1).m(); - } - - unsigned int n() const - { - return system_matrix->block(1,1).n(); - } - - private: - const SmartPointer > system_matrix; - const SmartPointer, Preconditioner> > A_inverse; - - mutable Vector tmp1, tmp2; - }; - - - - template - SchurComplement:: - SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix,Preconditioner> &A_inverse) - : - system_matrix (&system_matrix), - A_inverse (&A_inverse), - tmp1 (system_matrix.block(0,0).m()), - tmp2 (system_matrix.block(0,0).m()) - {} - - - template - void SchurComplement::vmult (Vector &dst, - const Vector &src) const + + template + void Obstacle::vector_value (const Point &p, + Vector &values) const { - system_matrix->block(0,1).vmult (tmp1, src); - A_inverse->name = "in schur"; - A_inverse->vmult (tmp2, tmp1); - system_matrix->block(1,0).vmult (dst, tmp2); - dst *= -1; - system_matrix->block(1,1).vmult_add (dst, src); - dst *= -1; + for (unsigned int c=0; cn_components; ++c) + values(c) = Obstacle::value (p, c); } +} + // @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 () + : + 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 (2000), + e_modul (2.e5), + 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)); +} - template - StokesProblem::StokesProblem (const unsigned int degree) - : - degree (degree), - triangulation (Triangulation::limit_level_difference_at_vertices), - fe (FE_Q(degree+1), dim, - FE_Q(degree), 1), - dof_handler (triangulation) - {} - - - - +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| + | 7 | / + |_______|/ + 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)->set_boundary_indicator (8); + if (cell->face (face)->center ()[1] == p1(1) || + cell->face (face)->center ()[1] == p2(1)) + cell->face (face)->set_boundary_indicator (7); + if (cell->face (face)->center ()[2] == p1(2)) + cell->face (face)->set_boundary_indicator (6); + } + + n_refinements = 2; + n_refinements_local = 3; + triangulation.refine_global (n_refinements); + + // 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) + // Verfeinerung fuer Dortmund + // if(// Einschraenkung in x-Richtung linker Rand + // cell->face (face)->vertex (0)[0] <= 0.7 &&//p2(0)*0.5+0.4 - 0.4*hlp_refinement && + // // Einschraenkung in x-Richtung rechter Rand + // cell->face (face)->vertex (1)[0] >= 0.3 &&//p2(0)*0.5-0.3 + 0.2*hlp_refinement && + // // Einschraenkung in y-Richtung oberer Rand + // cell->face (face)->vertex (0)[1] <= 0.875 &&//p2(1)-0.6 + 0.4*hlp_refinement) + // // Einschraenkung in y-Richtung unterer Rand + // cell->face (face)->vertex (2)[1] >= 0.125) + // { + // cell->set_refine_flag (); + // break; + // } + // Verfeinerung TKSE + { + cell->set_refine_flag (); + break; + } + else if (cell->level () == n_refinements + n_refinements_local - 1) + { + cell->set_refine_flag (); + break; + } + }; + + // Markierte Zellen werden verfeinert + triangulation.execute_coarsening_and_refinement (); + }; +} - template - void StokesProblem::setup_dofs () +template +void Step4::setup_system () +{ + // setup dofs { - A_preconditioner.reset (); - mg_A_preconditioner.resize (0); - system_matrix.clear (); - dof_handler.distribute_dofs (fe); -// DoFRenumbering::Cuthill_McKee (dof_handler); - std::vector block_component (dim+1,0); - block_component[dim] = 1; - DoFRenumbering::component_wise (dof_handler, block_component); - - - { - constraints.clear (); - typename FunctionMap::type dirichlet_boundary; - ZeroFunction homogeneous_dirichlet_bc (dim+1); //TODO: go back to BoundaryValues - - dirichlet_boundary[0] = &homogeneous_dirichlet_bc; - MappingQ1 mapping; - - std::vector component_mask (dim+1, true); - component_mask[dim] = false; - VectorTools::interpolate_boundary_values (mapping, - dof_handler, - dirichlet_boundary, - constraints, - component_mask); - - DoFTools::make_hanging_node_constraints (dof_handler, - constraints); - - mg_constrained_dofs.clear(); - mg_constrained_dofs.initialize(dof_handler, dirichlet_boundary); - } - - constraints.close (); - - - std::vector dofs_per_block (2); - DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, - block_component); - const unsigned int n_u = dofs_per_block[0], - n_p = dofs_per_block[1]; - - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl - << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << " (" << n_u << '+' << n_p << ')' - << std::endl; - - { - BlockCompressedSimpleSparsityPattern csp (2,2); - - csp.block(0,0).reinit (n_u, n_u); - csp.block(1,0).reinit (n_p, n_u); - csp.block(0,1).reinit (n_u, n_p); - csp.block(1,1).reinit (n_p, n_p); - - csp.collect_sizes(); - - DoFTools::make_sparsity_pattern ( - static_cast&>(dof_handler), - csp, constraints, false); - sparsity_pattern.copy_from (csp); - } - - - system_matrix.reinit (sparsity_pattern); - - solution.reinit (2); - solution.block(0).reinit (n_u); - solution.block(1).reinit (n_p); - solution.collect_sizes (); - - system_rhs.reinit (2); - system_rhs.block(0).reinit (n_u); - system_rhs.block(1).reinit (n_p); - system_rhs.collect_sizes (); - - //now setup stuff for mg - const unsigned int nlevels = triangulation.n_levels(); - - mg_matrices.resize(0, nlevels-1); - mg_matrices.clear (); - mg_interface_matrices.resize(0, nlevels-1); - mg_interface_matrices.clear (); - mg_sparsity.resize(0, nlevels-1); - - mg_dofs_per_component.resize (nlevels); - for (unsigned int level=0; level - void StokesProblem::assemble_system () + // setup hanging nodes and dirichlet constraints { - system_matrix=0; - system_rhs=0; - - QGauss quadrature_formula(degree+2); - - FEValues fe_values (fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - update_gradients); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - Vector local_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const RightHandSide right_hand_side; - std::vector > rhs_values (n_q_points, - Vector(dim+1)); - - - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); - - - - std::vector > phi_grads_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); - - typename MGDoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - local_matrix = 0; - local_rhs = 0; - - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - rhs_values); - - for (unsigned int q=0; qget_dof_indices (local_dof_indices); - constraints.distribute_local_to_global (local_matrix, local_rhs, - local_dof_indices, - system_matrix, system_rhs); - } + // 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 (); } - - template - void StokesProblem::assemble_multigrid () + // Initialzation for matrices and vectors { - QGauss quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - update_gradients); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); - - - std::vector > phi_grads_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); - - std::vector > interface_dofs - = mg_constrained_dofs.get_refinement_edge_indices (); - std::vector > boundary_interface_dofs - = mg_constrained_dofs.get_refinement_edge_boundary_indices (); - - std::vector boundary_constraints (triangulation.n_levels()); - std::vector boundary_interface_constraints (triangulation.n_levels()); - for (unsigned int level=0; level::cell_iterator - cell = dof_handler.begin(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - // Remember the level of the - // current cell. - const unsigned int level = cell->level(); - // Compute the values specified - // by update flags above. - fe_values.reinit (cell); - local_matrix = 0; - - for (unsigned int q=0; qget_mg_dof_indices (local_dof_indices); - boundary_constraints[level] - .distribute_local_to_global (local_matrix, - local_dof_indices, - mg_matrices[level]); - - for (unsigned int i=0; i::type>(new typename InnerPreconditioner::type()); - mg_A_preconditioner[level] - ->initialize (mg_matrices[level].block(0,0), - typename InnerPreconditioner::type::AdditionalData()); - } + 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 ()); } - - template - class SchurComplementSmoother + // setup sparsity pattern { - public: - struct AdditionalData - { - const InnerPreconditioner *A_preconditioner; - }; - - void initialize (const BlockSparseMatrix &system_matrix, - const AdditionalData &data); - - void vmult (BlockVector &dst, - const BlockVector &src) const; - - void Tvmult (BlockVector &dst, - const BlockVector &src) const; - - void clear (); - - private: - SmartPointer > system_matrix; - SmartPointer A_preconditioner; - }; - - - template - void - SchurComplementSmoother:: - initialize (const BlockSparseMatrix &system_matrix, - const AdditionalData &data) - { - this->system_matrix = &system_matrix; - this->A_preconditioner = data.A_preconditioner; + 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); + + mass_matrix.reinit (sp); } - - - - template - void - SchurComplementSmoother:: - vmult (BlockVector &dst, - const BlockVector &src) const - { -#ifdef STEP_42_TEST - std::cout << "Entering smoother with " << dst.size() << " unknowns" << std::endl; -#endif - - SparseDirectUMFPACK direct_solver; - direct_solver.initialize(*system_matrix); - Vector solution, rhs; - solution = dst; - rhs = src; - direct_solver.vmult(solution, rhs); - dst = solution; -/* - const InverseMatrix,InnerPreconditioner> - A_inverse (system_matrix->block(0,0), *A_preconditioner); - Vector tmp (dst.block(0).size()); - - - { - Vector schur_rhs (dst.block(1).size()); - A_inverse.name = "rhs"; - A_inverse.vmult (tmp, src.block(0)); -// std::cout << " TMP " << tmp.l2_norm() << std::endl; -system_matrix->block(1,0).vmult (schur_rhs, tmp); -schur_rhs -= src.block(1); -// std::cout << " BLOCK 1 " << src.block(1).l2_norm() << std::endl; -// std::cout << " SCHUR RHS " << schur_rhs.l2_norm() << std::endl; - -SchurComplement -schur_complement (*system_matrix, A_inverse); - - // The usual control structures for - // the solver call are created... - SolverControl solver_control (dst.block(1).size(), - 1e-1*schur_rhs.l2_norm()); - SolverGMRES<> cg (solver_control); - - #ifdef STEP_42_TEST - std::cout << " Starting Schur complement solver -- " - << schur_complement.m() << " unknowns" - << std::endl; - #endif - try - { - cg.solve (schur_complement, dst.block(1), schur_rhs, - PreconditionIdentity()); - } - catch (...) - { - std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; - std::cout << schur_rhs.l2_norm () << std::endl; - abort (); - } - -// no constraints to be taken care of here -#ifdef STEP_42_TEST -std::cout << " " -<< solver_control.last_step() -<< " CG Schur complement iterations in smoother " -<< solver_control.initial_value() << "->" << solver_control.last_value() -<< std::endl; -#endif + assemble_mass_matrix (); + const unsigned int + start = (system_rhs_newton.local_range().first), + end = (system_rhs_newton.local_range().second); + for (unsigned int j=0; j +void Step4::assemble_mass_matrix () +{ + QTrapez face_quadrature_formula; + + 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 dofs_per_face = fe.dofs_per_face; + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, 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(); + + 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 (); +} +template +void Step4::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) { -system_matrix->block(0,1).vmult (tmp, dst.block(1)); -tmp *= -1; -tmp += src.block(0); + QGauss quadrature_formula(2); + QGauss face_quadrature_formula(2); -A_inverse.name = "top left"; -A_inverse.vmult (dst.block(0), tmp); -// no constraints here either -} -#ifdef STEP_42_TEST -std::cout << "Exiting smoother with " << dst.size() << " unknowns" << std::endl; -#endif -*/ - } + 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(); - template - void - SchurComplementSmoother::clear () - {} + 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); - template - void - SchurComplementSmoother:: - Tvmult (BlockVector &, - const BlockVector &) const - { - Assert (false, ExcNotImplemented()); - } - + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + const FEValuesExtractors::Vector displacement (0); - template - void StokesProblem::solve () - { - system_matrix.block(1,1) = 0; - assemble_multigrid (); - typedef PreconditionMG, MGTransferPrebuilt > > - MGPREC; - - GrowingVectorMemory > mg_vector_memory; - - MGTransferPrebuilt > mg_transfer(constraints, mg_constrained_dofs); - std::vector block_component (dim+1,0); - block_component[dim] = 1; - mg_transfer.set_component_to_block_map (block_component); - mg_transfer.build_matrices(dof_handler); - - FullMatrix mg_coarse_matrix; - mg_coarse_matrix.copy_from (mg_matrices[0]); - MGCoarseGridHouseholder > mg_coarse; - mg_coarse.initialize(mg_coarse_matrix); - - MGMatrix, BlockVector > - mg_matrix(&mg_matrices); - MGMatrix, BlockVector > - mg_interface_up(&mg_interface_matrices); - MGMatrix, BlockVector > - mg_interface_down(&mg_interface_matrices); - - typedef - SchurComplementSmoother::type> - Smoother; - - MGSmootherPrecondition, - Smoother, - BlockVector > - mg_smoother(mg_vector_memory); - - MGLevelObject - smoother_data (0, triangulation.n_levels()-1); - - for (unsigned int level=0; level > mg(dof_handler, - mg_matrix, - mg_coarse, - mg_transfer, - mg_smoother, - mg_smoother); - mg.set_debug(3); - mg.set_edge_matrices(mg_interface_down, mg_interface_up); - - MGPREC preconditioner(dof_handler, mg, mg_transfer); - - SolverControl solver_control (system_matrix.m(), - 1e-6*system_rhs.l2_norm()); - GrowingVectorMemory > vector_memory; - SolverGMRES >::AdditionalData gmres_data; - gmres_data.max_n_tmp_vectors = 100; - - SolverGMRES > gmres(solver_control, vector_memory, - gmres_data); - -// PreconditionIdentity precondition_identity; -#ifdef STEP_42_TEST - std::cout << "Starting outer GMRES complement solver" << std::endl; -#endif - try - { - gmres.solve(system_matrix, solution, system_rhs, - preconditioner); - } - catch (...) + TrilinosWrappers::MPI::Vector test_rhs(solution); + const double kappa = 1.0; + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) { - std::cout << "Failure in " << __PRETTY_FUNCTION__ << std::endl; - abort (); - } - - constraints.distribute (solution); - - std::cout << solver_control.last_step() - << " outer GMRES iterations "; - } - - - template - void StokesProblem::solve_block () - { - std::cout << " Computing preconditioner..." << std::endl << std::flush; - - A_preconditioner - = std_cxx1x::shared_ptr::type>(new typename InnerPreconditioner::type()); - A_preconditioner->initialize (system_matrix.block(0,0), - typename InnerPreconditioner::type::AdditionalData()); - - SparseMatrix pressure_mass_matrix; - pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); - pressure_mass_matrix.copy_from(system_matrix.block(1,1)); - system_matrix.block(1,1) = 0; - - SparseILU pmass_preconditioner; - pmass_preconditioner.initialize (pressure_mass_matrix, - SparseILU::AdditionalData()); - - InverseMatrix,SparseILU > - m_inverse (pressure_mass_matrix, pmass_preconditioner); - - BlockSchurPreconditioner::type, - SparseILU > - preconditioner (system_matrix, m_inverse, *A_preconditioner); - - SolverControl solver_control (system_matrix.m(), - 1e-6*system_rhs.l2_norm()); - GrowingVectorMemory > vector_memory; - SolverGMRES >::AdditionalData gmres_data; - gmres_data.max_n_tmp_vectors = 100; - - SolverGMRES > gmres(solver_control, vector_memory, - gmres_data); - - gmres.solve(system_matrix, solution, system_rhs, - preconditioner); - - constraints.distribute (solution); + 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;//right_hand_side_values[q_point](0); + //rhs_values[1] = 0;//right_hand_side_values[q_point](1); + cell_rhs(i) += (fe_values[displacement].value (i, q_point) * + rhs_values * + fe_values.JxW (q_point)); + + // if (i == 7) + // std::cout<< i << ". " << stress_tensor//->get_strain(fe_values, i, q_point) + // << ", " << stress_strain_tensor_linearized + // <::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] = right_hand_side_values_face[q_point](0); + rhs_values[1] = right_hand_side_values_face[q_point](1); + rhs_values[2] = right_hand_side_values_face[q_point](2); + + 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 (); + + // for (unsigned int i=0; i +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] = right_hand_side_values[q_point](0); + rhs_values[1] = right_hand_side_values[q_point](1); + 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] = right_hand_side_values_face[q_point](0); + rhs_values[1] = right_hand_side_values_face[q_point](1); + rhs_values[2] = right_hand_side_values_face[q_point](2); + + 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); + // for (unsigned int i=0; i - void - StokesProblem::output_results (const unsigned int refinement_cycle) const - { - std::vector solution_names (dim, "velocity"); - solution_names.push_back ("pressure"); - - std::vector - data_component_interpretation - (dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - - DataOut data_out; - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, solution_names, - DataOut::type_dof_data, - data_component_interpretation); - data_out.build_patches (); - - std::ostringstream filename; - filename << "solution-" - << Utilities::int_to_string (refinement_cycle, 2) - << ".vtk"; - - std::ofstream output (filename.str().c_str()); - data_out.write_vtk (output); - } + 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 < 0) + { + 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.add_index (index_z); + + // std::cout<< point[0] << " " + // << -lambda (index_z) << " " + // << diag_mass_matrix_vector (index_z) << " " + // << 1.0/pow (2.0,n_refinements + n_refinements_local) + // < - void - StokesProblem::refine_mesh () - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - std::vector component_mask (dim+1, false); - component_mask[dim] = true; - KellyErrorEstimator::estimate (static_cast&>(dof_handler), - QGauss(degree+1), - typename FunctionMap::type(), - solution, - estimated_error_per_cell, - component_mask); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.0); - triangulation.execute_coarsening_and_refinement (); - } +template +void Step4::dirichlet_constraints () +{ + /* boundary_indicators: + _______ + / 9 /| + /______ / | + 8| | 8| + | 7 | / + |_______|/ + 6 + */ + + // constraints_dirichlet_hanging_nodes.clear (); + 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] = false; + 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] = false; + component_mask[1] = true; + component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 7, + EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, + component_mask); + + component_mask[0] = true; + component_mask[1] = false; + component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 8, + EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, + component_mask); + constraints_dirichlet_hanging_nodes.close (); +} +template +void Step4::solve () +{ + ReductionControl reduction_control (10000, 1e-15, 1e-4); + + TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); + distributed_solution = solution; + + // constraints_hanging_nodes.set_zero (distributed_solution); + const unsigned int + start = (distributed_solution.local_range().first), + end = (distributed_solution.local_range().second); + for (unsigned int i=start; i + solver (reduction_control, mpi_communicator); + + preconditioner_u.initialize (system_matrix_newton, additional_data); + // preconditioner_u.reinit (); + + solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); + pcout << "Initial error: " << reduction_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); + clock_t start, end; + + char* name = new char[30]; + sprintf (name,"Error-FPV.dat"); + FILE* fp = fopen(name,"w"); + delete[] name; + fclose(fp); + + 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; + //additional_data.constant_modes; + + // std::cout<< "Update Active Set in Dim = " << dim < 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<< " " < " << resid < subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i subdivisions (dim, 1); - subdivisions[0] = 1; - - const Point bottom_left = (dim == 2 ? - Point(0,0) : - Point(0,0,0)); - const Point top_right = (dim == 2 ? - Point(1,1) : - Point(1,1,1)); - - GridGenerator::subdivided_hyper_rectangle (triangulation, - subdivisions, - bottom_left, - top_right); + std::vector filenames; + for (unsigned int i=0; + i +void Step4::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; + } + } +} - for (typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->center()[0] == 1) - cell->face(f)->set_all_boundary_indicators(1); - - - - triangulation.refine_global (1); +template +void Step4::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"); - for (unsigned int refinement_cycle = 0; refinement_cycle<10; - ++refinement_cycle) - { - std::cout << "Refinement cycle " << refinement_cycle << std::endl; + data_out.build_patches (); - if (refinement_cycle > 0) - refine_mesh (); + std::ofstream output_vtk (dim == 2 ? + (title + ".vtk").c_str () : + (title + ".vtk").c_str ()); + data_out.write_vtk (output_vtk); +} - std::ostringstream out_filename; - out_filename << "gitter" - << refinement_cycle - << ".eps"; +template +void Step4::output_results (Vector vector, const std::string& title) const +{ + DataOut data_out; - std::ofstream grid_output (out_filename.str().c_str()); - GridOut grid_out; - grid_out.write_eps (triangulation, grid_output); + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (vector, "vector_to_plot"); - setup_dofs (); + data_out.build_patches (); - std::cout << " Assembling..." << std::endl << std::flush; - assemble_system (); + std::ofstream output_vtk (dim == 2 ? + (title + ".vtk").c_str () : + (title + ".vtk").c_str ()); + data_out.write_vtk (output_vtk); +} - std::cout << " Solving..." << std::flush; +template +void Step4::run () +{ + pcout << "Solving problem in " << dim << " space dimensions." << std::endl; - solve_block (); - output_results (refinement_cycle); - system ("mv solution-* block"); + run_time.resize (5); - solution = 0; + clock_t start, end; - solve (); - output_results (refinement_cycle); - system ("mv solution-* mg"); + start = clock(); + make_grid(); + // mesh_surface (); + setup_system (); + end = clock(); + run_time[0] = (double)(end-start)/CLOCKS_PER_SEC; - std::cout << std::endl; - } - } + solve_newton (); } - -int main () + // @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. + // + // Finally, the first line of the function is + // used to suppress some output. Remember + // that in the previous example, we had the + // output from the linear solvers about the + // starting residual and the number of the + // iteration where convergence was + // detected. This can be suppressed through + // the deallog.depth_console(0) + // call. + // + // The rationale here is the following: the + // deallog (i.e. deal-log, not de-allog) + // variable represents a stream to which some + // parts of the library write output. It + // redirects this output to the console and + // if required to a file. The output is + // nested in a way so that each function can + // use a prefix string (separated by colons) + // for each line of output; if it calls + // another function, that may also use its + // prefix which is then printed after the one + // of the calling function. Since output from + // functions which are nested deep below is + // usually not as important as top-level + // output, you can give the deallog variable + // a maximal depth of nested output for + // output to console and file. The depth zero + // which we gave here means that no output is + // written. By changing it you can get more + // information about the innards of the + // library. +int main (int argc, char *argv[]) { - try - { - using namespace dealii; - using namespace Step42; + deallog.depth_console (0); - deallog.depth_console (0); + clock_t start, end; - StokesProblem<2> flow_problem(1); - flow_problem.run (); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } + start = clock(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv); + { + Step4<3> laplace_problem_3d; + laplace_problem_3d.run (); + } + end = clock(); + cout<< "%%%%%% Rechenzeit overall = " << (double)(end-start)/CLOCKS_PER_SEC <