]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Initial version of elasto-plastic contact problem
authorfrohne <frohne@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Apr 2012 13:00:37 +0000 (13:00 +0000)
committerfrohne <frohne@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Apr 2012 13:00:37 +0000 (13:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@25356 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-42/Makefile
deal.II/examples/step-42/step-42.cc

index 8ffa568be9871fcae131aefc0a031f7f872cd4d2..4f520f02a301c656b5c9dae2c62e864402e04c48 100644 (file)
@@ -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
index 0c5fab41704e83c9b7d988f44d6f33269b7ee02b..c4e9230ad1e536c658deb08196399fad91de115a 100644 (file)
-/* 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 <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/constraint_matrix.h>
+                                // The first few (many?) include
+                                // files have already been used in
+                                // the previous example, so we will
+                                // not explain their meaning here
+                                // again.
 
 #include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/grid_refinement.h>
-
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_renumbering.h>
 #include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-
+#include <deal.II/dofs/dof_renumbering.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_system.h>
+#include <deal.II/dofs/dof_tools.h>
 #include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping_q1.h>
-
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
 #include <deal.II/numerics/vectors.h>
 #include <deal.II/numerics/matrices.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/error_estimator.h>
-
-
-#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
 
-#include <deal.II/lac/sparse_ilu.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_block_vector.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
 
-#include <deal.II/multigrid/multigrid.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
-#include <deal.II/multigrid/mg_dof_accessor.h>
-#include <deal.II/multigrid/mg_constrained_dofs.h>
-#include <deal.II/multigrid/mg_transfer.h>
-#include <deal.II/multigrid/mg_tools.h>
-#include <deal.II/multigrid/mg_coarse.h>
-#include <deal.II/multigrid/mg_smoother.h>
-#include <deal.II/multigrid/mg_matrix.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
 
+#include <deal.II/numerics/data_out.h>
 #include <fstream>
-#include <sstream>
-
-
-namespace Step42
-{
-  using namespace dealii;
+#include <iostream>
+#include <list>
+#include <time.h>
+
+                                // 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 <deal.II/base/logstream.h>
 
+                                // 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 <int dim>
-  struct InnerPreconditioner;
+                                 // @sect3{The <code>Step4</code> class template}
 
+template <int dim> class ConstitutiveLaw;
 
-  template <>
-  struct InnerPreconditioner<2>
-  {
-      typedef SparseDirectUMFPACK type;
-  };
+template <int dim>
+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<double>                &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<double> vector, const std::string& title) const;
+
+  MPI_Comm             mpi_communicator;  
+
+  parallel::distributed::Triangulation<dim>   triangulation;
+
+  FESystem<dim>        fe;
+  DoFHandler<dim>      dof_handler;
+  
+  IndexSet             locally_owned_dofs;
+  IndexSet             locally_relevant_dofs;
+  
+  int                  n_refinements;
+  int                  n_refinements_local;
+  unsigned int         number_iterations;
+  std::vector<double>  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<ConstitutiveLaw<dim> > 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<TrilinosWrappers::PreconditionAMG>  Mp_preconditioner;
+};
+
+template <int dim>
+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<dim> &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 <int dim>
+ConstitutiveLaw<dim>::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 <<std::endl;
+  stress_strain_tensor_kappa = kappa*outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>());
+  stress_strain_tensor_mu = 2*mu*(identity_tensor<dim>() - outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>())/3.0);
+}
 
+template <int dim>
+inline
+SymmetricTensor<2,dim> ConstitutiveLaw<dim>::get_strain (const FEValues<dim> &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<double> type;
-  };
+  tmp = fe_values[displacement].symmetric_gradient (shape_func,q_point);
+  
+  return tmp;
+}
 
-  template <typename MATRIX>
-  void copy(const MATRIX &matrix,
-            FullMatrix<double> &full_matrix)
+template <int dim>
+void ConstitutiveLaw<dim>::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 <<std::endl;
+    stress_strain_tensor += stress_strain_tensor_kappa;
 
-    Vector<double> unit (n);
-    Vector<double> result (m);
-    for(unsigned int i=0; i<n; ++i)
-      {
-        unit(i) = 1;
-        for(unsigned int j=0; j<m; ++j)
-          {
-            matrix.vmult(result,unit);
-            full_matrix(i,j) = result(j);
-          }
-        unit(i) = 0;
-      }
+    sigma_eff = beta * deviator_stress_tensor_norm;
   }
+}
 
-  template <int dim>
-  class StokesProblem
+template <int dim>
+void ConstitutiveLaw<dim>::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<std::vector<bool> > &lower_dofs,
-                                     std::vector<std::vector<bool> > &boundary_dofs);
-
-      void output_results (const unsigned int refinement_cycle) const;
-      void refine_mesh ();
-
-      const unsigned int   degree;
-
-      Triangulation<dim>   triangulation;
-      FESystem<dim>        fe;
-      MGDoFHandler<dim>    dof_handler;
-
-      ConstraintMatrix     constraints;
-
-      BlockSparsityPattern      sparsity_pattern;
-      BlockSparseMatrix<double> system_matrix;
-
-      BlockVector<double> solution;
-      BlockVector<double> system_rhs;
-
-      MGLevelObject<ConstraintMatrix>           mg_constraints;
-      MGLevelObject<BlockSparsityPattern>       mg_sparsity;
-      MGLevelObject<BlockSparseMatrix<double> > mg_matrices;
-
-      MGLevelObject<BlockSparseMatrix<double> > mg_interface_matrices;
-      MGConstrainedDoFs                         mg_constrained_dofs;
-      std::vector<std::vector<unsigned int> >   mg_dofs_per_component;
-
-      std::vector<std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type> > mg_A_preconditioner;
-      std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::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 <int dim>
-  class BoundaryValues : public Function<dim>
+  class RightHandSide : public Function<dim> 
   {
-    public:
-      BoundaryValues () : Function<dim>(dim+1) {}
-
-      virtual double value (const Point<dim>   &p,
-                            const unsigned int  component = 0) const;
-
-      virtual void vector_value (const Point<dim> &p,
-                                 Vector<double>   &value) const;
+  public:
+    RightHandSide () : Function<dim>(dim) {}
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+    
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double>   &values) const;
   };
 
-
   template <int dim>
-  double
-  BoundaryValues<dim>::value (const Point<dim>  &p,
-                              const unsigned int component) const
+  double RightHandSide<dim>::value (const Point<dim> &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<dim; ++i)
+    //   return_value += 4*std::pow(p(i), 4);
+    
+    return return_value;
   }
-
-
+  
   template <int dim>
-  void
-  BoundaryValues<dim>::vector_value (const Point<dim> &p,
-                                     Vector<double>   &values) const
+  void RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                        Vector<double>   &values) const
   {
     for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = BoundaryValues<dim>::value (p, c);
+      values(c) = RightHandSide<dim>::value (p, c);
   }
 
 
-
-
   template <int dim>
-  class RightHandSide : public Function<dim>
+  class BoundaryValues : public Function<dim> 
   {
-    public:
-      RightHandSide () : Function<dim>(dim+1) {}
-
-      virtual double value (const Point<dim>   &p,
-                            const unsigned int  component = 0) const;
-
-      virtual void vector_value (const Point<dim> &p,
-                                 Vector<double>   &value) const;
-
+  public:
+    BoundaryValues () : Function<dim>(dim) {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+    
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double>   &values) const;
   };
-
-
+  
   template <int dim>
-  double
-  RightHandSide<dim>::value (const Point<dim>  &/*p*/,
-                             const unsigned int component) const
+  double BoundaryValues<dim>::value (const Point<dim> &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 <int dim>
-  void
-  RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                    Vector<double>   &values) const
+  void BoundaryValues<dim>::vector_value (const Point<dim> &p,
+                                          Vector<double>   &values) const 
   {
     for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = RightHandSide<dim>::value (p, c);
+      values(c) = BoundaryValues<dim>::value (p, c);
   }
 
-
-
-
-  template <class Matrix, class Preconditioner>
-  class InverseMatrix : public Subscriptor
+  
+  template <int dim>
+  class Obstacle : public Function<dim> 
   {
-    public:
-      InverseMatrix (const Matrix         &m,
-                     const Preconditioner &preconditioner);
+  public:
+    Obstacle () : Function<dim>(dim) {};
 
-      void vmult (Vector<double>       &dst,
-                  const Vector<double> &src) const;
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
 
-      mutable std::string name;
-    private:
-      const SmartPointer<const Matrix> matrix;
-      const SmartPointer<const Preconditioner> preconditioner;
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double>   &values) const;
   };
 
-
-  template <class Matrix, class Preconditioner>
-  InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
-                                                       const Preconditioner &preconditioner)
-                  :
-                  matrix (&m),
-                  preconditioner (&preconditioner)
-  {}
-
-
-  template <class Matrix, class Preconditioner>
-  void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
-                                                    const Vector<double> &src) const
+  template <int dim>
+  double Obstacle<dim>::value (const Point<dim> &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 PreconditionerA, class PreconditionerMp>
-  class BlockSchurPreconditioner : public Subscriptor
-  {
-    public:
-      BlockSchurPreconditioner (const BlockSparseMatrix<double>         &S,
-                                const InverseMatrix<SparseMatrix<double>,PreconditionerMp>  &Mpinv,
-                                const PreconditionerA &Apreconditioner);
-
-      void vmult (BlockVector<double>       &dst,
-                  const BlockVector<double> &src) const;
-
-    private:
-      const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-      const SmartPointer<const InverseMatrix<SparseMatrix<double>,
-                                             PreconditionerMp > > m_inverse;
-      const PreconditionerA &a_preconditioner;
-
-      mutable Vector<double> tmp;
-
-  };
-
-  template <class PreconditionerA, class PreconditionerMp>
-  BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::BlockSchurPreconditioner(
-    const BlockSparseMatrix<double>                            &S,
-    const InverseMatrix<SparseMatrix<double>,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 <class PreconditionerA, class PreconditionerMp>
-  void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
-    BlockVector<double>       &dst,
-    const BlockVector<double> &src) const
-  {
-                                     // Form u_new = A^{-1} u
-    a_preconditioner.vmult (dst.block(0), src.block(0));
-                                     // Form tmp = - B u_new + p
-                                     // (<code>SparseMatrix::residual</code>
-                                     // 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 Preconditioner>
-  class SchurComplement : public Subscriptor
-  {
-    public:
-      SchurComplement (const BlockSparseMatrix<double> &system_matrix,
-                       const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
-
-      void vmult (Vector<double>       &dst,
-                  const Vector<double> &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<const BlockSparseMatrix<double> > system_matrix;
-      const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
-
-      mutable Vector<double> tmp1, tmp2;
-  };
-
-
-
-  template <class Preconditioner>
-  SchurComplement<Preconditioner>::
-  SchurComplement (const BlockSparseMatrix<double> &system_matrix,
-                   const InverseMatrix<SparseMatrix<double>,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 <class Preconditioner>
-  void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
-                                               const Vector<double> &src) const
+  
+  template <int dim>
+  void Obstacle<dim>::vector_value (const Point<dim> &p,
+                                   Vector<double> &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; c<this->n_components; ++c)
+      values(c) = Obstacle<dim>::value (p, c);  
   }
+}
 
 
+                                 // @sect3{Implementation of the <code>Step4</code> class}
+
+                                 // Next for the implementation of the class
+                                 // template that makes use of the functions
+                                 // above. As before, we will write everything
+              
+template <int dim>
+Step4<dim>::Step4 ()
+  :
+  mpi_communicator (MPI_COMM_WORLD),
+  triangulation (mpi_communicator),
+  fe (FE_Q<dim>(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<dim> (e_modul, nu, sigma_0, gamma, mpi_communicator, pcout));
+}
 
-  template <int dim>
-  StokesProblem<dim>::StokesProblem (const unsigned int degree)
-                  :
-                  degree (degree),
-                  triangulation (Triangulation<dim>::limit_level_difference_at_vertices),
-                  fe (FE_Q<dim>(degree+1), dim,
-                      FE_Q<dim>(degree), 1),
-                  dof_handler (triangulation)
-  {}
-
-
-
-
+template <int dim>
+void Step4<dim>::make_grid ()
+{
+  std::vector<unsigned int> repet(3);
+  repet[0] = 1;//20;
+  repet[1] = 1;
+  repet[2] = 1;
+
+  Point<dim> p1 (0,0,0);
+  Point<dim> 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<GeometryInfo<dim>::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<n_refinements_local; ++step)
+    {
+      cell = triangulation.begin_active();  // Iterator ueber alle Zellen
+     
+      double hlp_refinement = 0;
+      hlp_refinement = pow((double)(step)/(n_refinements_local),4.0);
+      pcout<< "Verfeinerungsfaktor: " << hlp_refinement <<std::endl;
+
+      for (; cell!=endc; ++cell)
+         for (unsigned int face=0; face<GeometryInfo<dim>::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 <int dim>
-  void StokesProblem<dim>::setup_dofs ()
+template <int dim>
+void Step4<dim>::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<unsigned int> block_component (dim+1,0);
-    block_component[dim] = 1;
-    DoFRenumbering::component_wise (dof_handler, block_component);
-
-
-    {
-      constraints.clear ();
-      typename FunctionMap<dim>::type      dirichlet_boundary;
-      ZeroFunction<dim>                    homogeneous_dirichlet_bc (dim+1); //TODO: go back to BoundaryValues
-
-      dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
-      MappingQ1<dim> mapping;
-
-      std::vector<bool> 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<unsigned int> 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<const DoFHandler<dim>&>(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<nlevels; ++level)
-      mg_dofs_per_component[level].resize (2);
-
-    MGTools::count_dofs_per_block (dof_handler, mg_dofs_per_component,
-                                   block_component);
-    for (unsigned int level=0; level<nlevels; ++level)
-      std::cout << "                        Level " << level << ": "
-                << dof_handler.n_dofs (level) << " ("
-                << mg_dofs_per_component[level][0] << '+'
-                << mg_dofs_per_component[level][1] << ')'
-                << std::endl;
-
-    for (unsigned int level=0; level<nlevels; ++level)
-      {
-        DoFRenumbering::component_wise (dof_handler, level, block_component);
-
-        BlockCompressedSparsityPattern bcsp (mg_dofs_per_component[level],
-                                             mg_dofs_per_component[level]);
-        MGTools::make_sparsity_pattern(dof_handler, bcsp, level);
-        mg_sparsity[level].copy_from (bcsp);
-        mg_matrices[level].reinit (mg_sparsity[level]);
-        mg_interface_matrices[level].reinit (mg_sparsity[level]);
-      }
+    locally_owned_dofs = dof_handler.locally_owned_dofs ();
+    locally_relevant_dofs.clear();
+    DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                            locally_relevant_dofs);
   }
 
-
-
-  template <int dim>
-  void StokesProblem<dim>::assemble_system ()
+  // setup hanging nodes and dirichlet constraints
   {
-    system_matrix=0;
-    system_rhs=0;
-
-    QGauss<dim>   quadrature_formula(degree+2);
-
-    FEValues<dim> 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<double>   local_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>       local_rhs (dofs_per_cell);
-
-    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-    const RightHandSide<dim>          right_hand_side;
-    std::vector<Vector<double> >      rhs_values (n_q_points,
-                                                  Vector<double>(dim+1));
-
-
-    const FEValuesExtractors::Vector velocities (0);
-    const FEValuesExtractors::Scalar pressure (dim);
-
-
-
-    std::vector<Tensor<2,dim> >          phi_grads_u (dofs_per_cell);
-    std::vector<double>                  div_phi_u   (dofs_per_cell);
-    std::vector<double>                  phi_p       (dofs_per_cell);
-
-    typename MGDoFHandler<dim>::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; q<n_q_points; ++q)
-          {
-            for (unsigned int k=0; k<dofs_per_cell; ++k)
-              {
-                phi_grads_u[k] = fe_values[velocities].gradient (k, q);
-                div_phi_u[k]   = fe_values[velocities].divergence (k, q);
-                phi_p[k]       = fe_values[pressure].value (k, q);
-              }
-
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              {
-                for (unsigned int j=0; j<dofs_per_cell; ++j)
-                  {
-                    local_matrix(i,j) += (scalar_product(phi_grads_u[i], phi_grads_u[j])
-                                          - div_phi_u[i] * phi_p[j]
-                                          - phi_p[i] * div_phi_u[j]
-                                          - phi_p[i] * phi_p[j]
-                    )
-                                         * fe_values.JxW(q);
-                  }
-
-                const unsigned int component_i =
-                  fe.system_to_component_index(i).first;
-                local_rhs(i) += fe_values.shape_value(i,q) *
-                                rhs_values[q](component_i) *
-                                fe_values.JxW(q);
-              }
-          }
-
-
-
-
-        cell->get_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 <int dim>
-  void StokesProblem<dim>::assemble_multigrid ()
+  // Initialzation for matrices and vectors
   {
-    QGauss<dim>   quadrature_formula(degree+2);
-    FEValues<dim> 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<double>   local_matrix (dofs_per_cell, dofs_per_cell);
-
-    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-    const FEValuesExtractors::Vector velocities (0);
-    const FEValuesExtractors::Scalar pressure (dim);
-
-
-    std::vector<Tensor<2,dim> >          phi_grads_u (dofs_per_cell);
-    std::vector<double>                  div_phi_u   (dofs_per_cell);
-    std::vector<double>                  phi_p       (dofs_per_cell);
-
-    std::vector<std::vector<bool> > interface_dofs
-      = mg_constrained_dofs.get_refinement_edge_indices ();
-    std::vector<std::vector<bool> > boundary_interface_dofs
-      = mg_constrained_dofs.get_refinement_edge_boundary_indices ();
-
-    std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
-    std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
-    for (unsigned int level=0; level<triangulation.n_levels(); ++level)
-      {
-        boundary_constraints[level].add_lines (interface_dofs[level]);
-        boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices()[level]);
-        boundary_constraints[level].close ();
-
-        boundary_interface_constraints[level]
-          .add_lines (boundary_interface_dofs[level]);
-        boundary_interface_constraints[level].close ();
-      }
-
-    typename MGDoFHandler<dim>::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; q<n_q_points; ++q)
-          {
-            for (unsigned int k=0; k<dofs_per_cell; ++k)
-              {
-                phi_grads_u[k] = fe_values[velocities].gradient (k, q);
-                div_phi_u[k]   = fe_values[velocities].divergence (k, q);
-                phi_p[k]       = fe_values[pressure].value (k, q);
-              }
-
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                local_matrix(i,j) += (
-                  scalar_product(phi_grads_u[i], phi_grads_u[j])
-                  - div_phi_u[i] * phi_p[j]
-                  - phi_p[i] * div_phi_u[j]
-//                                        - phi_p[i] * phi_p[j]
-                )
-                                     * fe_values.JxW(q);
-          }
-
-        cell->get_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<dofs_per_cell; ++i)
-          for (unsigned int j=0; j<dofs_per_cell; ++j)
-            if( !(interface_dofs[level][local_dof_indices[i]]==true &&
-                  interface_dofs[level][local_dof_indices[j]]==false))
-              local_matrix(i,j) = 0;
-
-        boundary_interface_constraints[level]
-          .distribute_local_to_global (local_matrix,
-                                       local_dof_indices,
-                                       mg_interface_matrices[level]);
-      }
-
-    mg_A_preconditioner.resize (triangulation.n_levels());
-    for (unsigned int level=0; level<triangulation.n_levels(); ++level)
-      {
-        mg_A_preconditioner[level]
-          = std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
-        mg_A_preconditioner[level]
-          ->initialize (mg_matrices[level].block(0,0),
-                        typename InnerPreconditioner<dim>::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 <typename InnerPreconditioner>
-  class SchurComplementSmoother
+  // setup sparsity pattern
   {
-    public:
-      struct AdditionalData
-      {
-          const InnerPreconditioner *A_preconditioner;
-      };
-
-      void initialize (const BlockSparseMatrix<double> &system_matrix,
-                       const AdditionalData            &data);
-
-      void vmult (BlockVector<double> &dst,
-                  const BlockVector<double> &src) const;
-
-      void Tvmult (BlockVector<double> &dst,
-                   const BlockVector<double> &src) const;
-
-      void clear ();
-
-    private:
-      SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-      SmartPointer<const InnerPreconditioner>        A_preconditioner;
-  };
-
-
-  template <typename InnerPreconditioner>
-  void
-  SchurComplementSmoother<InnerPreconditioner>::
-  initialize (const BlockSparseMatrix<double> &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 <typename InnerPreconditioner>
-  void
-  SchurComplementSmoother<InnerPreconditioner>::
-  vmult (BlockVector<double> &dst,
-         const BlockVector<double> &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<double> solution, rhs;
-    solution = dst;
-    rhs = src;
-    direct_solver.vmult(solution, rhs);
-    dst = solution;
-/*
-  const InverseMatrix<SparseMatrix<double>,InnerPreconditioner>
-  A_inverse (system_matrix->block(0,0), *A_preconditioner);
-  Vector<double> tmp (dst.block(0).size());
-
-
-  {
-  Vector<double> 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<InnerPreconditioner>
-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<end; j++)
+    diag_mass_matrix_vector (j) = mass_matrix.diag_element (j);
+  number_iterations = 0;
 }
 
+template <int dim>
+void Step4<dim>::assemble_mass_matrix () 
+{  
+  QTrapez<dim-1>  face_quadrature_formula;
+
+  FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula, 
+                                     update_values   | update_quadrature_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell      = fe.dofs_per_cell;
+  const unsigned int   dofs_per_face      = fe.dofs_per_face;
+  const unsigned int   n_face_q_points    = face_quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  const FEValuesExtractors::Vector displacement (0);
+
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    if (cell->is_locally_owned())
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+       if (cell->face (face)->at_boundary()
+           && cell->face (face)->boundary_indicator () == 9)
+         {
+           fe_values_face.reinit (cell, face);
+           cell_matrix = 0;
+           
+           for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               cell_matrix(i,i) += (fe_values_face[displacement].value (i, q_point) *
+                                    fe_values_face[displacement].value (i, q_point) *
+                                    fe_values_face.JxW (q_point));
+           
+           cell->get_dof_indices (local_dof_indices);
+           
+           constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix,
+                                                                           local_dof_indices,
+                                                                           mass_matrix);
+         }
+
+  mass_matrix.compress ();
+}
 
+template <int dim>
+void Step4<dim>::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) 
 {
-system_matrix->block(0,1).vmult (tmp, dst.block(1));
-tmp *= -1;
-tmp += src.block(0);
+  QGauss<dim>  quadrature_formula(2);
+  QGauss<dim-1>  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<dim> fe_values (fe, quadrature_formula, 
+                          UpdateFlags(update_values    |
+                                      update_gradients |
+                                      update_q_points  |
+                                      update_JxW_values));
 
+  FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula, 
+                                   update_values   | update_quadrature_points |
+                                   update_JxW_values);
 
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size ();
+  const unsigned int   n_face_q_points = face_quadrature_formula.size();
 
-  template <typename InnerPreconditioner>
-  void
-  SchurComplementSmoother<InnerPreconditioner>::clear ()
-  {}
+  const EquationData::RightHandSide<dim> right_hand_side;
+  std::vector<Vector<double> > right_hand_side_values (n_q_points,
+                                                       Vector<double>(dim));
+  std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
+                                                           Vector<double>(dim));
 
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);  
+  Vector<double>       cell_rhs (dofs_per_cell);
 
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-  template <typename InnerPreconditioner>
-  void
-  SchurComplementSmoother<InnerPreconditioner>::
-  Tvmult (BlockVector<double> &,
-          const BlockVector<double> &) const
-  {
-    Assert (false, ExcNotImplemented());
-  }
-
+  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                endc = dof_handler.end();
 
+  const FEValuesExtractors::Vector displacement (0);
 
-  template <int dim>
-  void StokesProblem<dim>::solve ()
-  {
-    system_matrix.block(1,1) = 0;
-    assemble_multigrid ();
-    typedef PreconditionMG<dim, BlockVector<double>, MGTransferPrebuilt<BlockVector<double> > >
-      MGPREC;
-
-    GrowingVectorMemory<BlockVector<double> >  mg_vector_memory;
-
-    MGTransferPrebuilt<BlockVector<double> > mg_transfer(constraints, mg_constrained_dofs);
-    std::vector<unsigned int> 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<float> mg_coarse_matrix;
-    mg_coarse_matrix.copy_from (mg_matrices[0]);
-    MGCoarseGridHouseholder<float, BlockVector<double> > mg_coarse;
-    mg_coarse.initialize(mg_coarse_matrix);
-
-    MGMatrix<BlockSparseMatrix<double>, BlockVector<double> >
-      mg_matrix(&mg_matrices);
-    MGMatrix<BlockSparseMatrix<double>, BlockVector<double> >
-      mg_interface_up(&mg_interface_matrices);
-    MGMatrix<BlockSparseMatrix<double>, BlockVector<double> >
-      mg_interface_down(&mg_interface_matrices);
-
-    typedef
-      SchurComplementSmoother<typename InnerPreconditioner<dim>::type>
-      Smoother;
-
-    MGSmootherPrecondition<BlockSparseMatrix<double>,
-      Smoother,
-      BlockVector<double> >
-      mg_smoother(mg_vector_memory);
-
-    MGLevelObject<typename Smoother::AdditionalData>
-      smoother_data (0, triangulation.n_levels()-1);
-
-    for (unsigned int level=0; level<triangulation.n_levels(); ++level)
-      smoother_data[level].A_preconditioner = mg_A_preconditioner[level].get();
-
-    mg_smoother.initialize(mg_matrices, smoother_data);
-    mg_smoother.set_steps(2);
-
-    Multigrid<BlockVector<double> > 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<BlockVector<double> > vector_memory;
-    SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
-    gmres_data.max_n_tmp_vectors = 100;
-
-    SolverGMRES<BlockVector<double> > 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 <int dim>
-  void StokesProblem<dim>::solve_block ()
-  {
-    std::cout << "   Computing preconditioner..." << std::endl << std::flush;
-
-    A_preconditioner
-      = std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
-    A_preconditioner->initialize (system_matrix.block(0,0),
-                                  typename InnerPreconditioner<dim>::type::AdditionalData());
-
-    SparseMatrix<double> 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<double> pmass_preconditioner;
-    pmass_preconditioner.initialize (pressure_mass_matrix,
-                                     SparseILU<double>::AdditionalData());
-
-    InverseMatrix<SparseMatrix<double>,SparseILU<double> >
-      m_inverse (pressure_mass_matrix, pmass_preconditioner);
-
-    BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
-      SparseILU<double> >
-      preconditioner (system_matrix, m_inverse, *A_preconditioner);
-
-    SolverControl solver_control (system_matrix.m(),
-                                  1e-6*system_rhs.l2_norm());
-    GrowingVectorMemory<BlockVector<double> > vector_memory;
-    SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
-    gmres_data.max_n_tmp_vectors = 100;
-
-    SolverGMRES<BlockVector<double> > 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<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
+       fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
+       
+       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+         {     
+           SymmetricTensor<4,dim> stress_strain_tensor_linearized;
+           SymmetricTensor<4,dim> stress_strain_tensor;
+           SymmetricTensor<2,dim> stress_tensor;
+           
+           plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized,
+                                                              stress_strain_tensor,
+                                                              strain_tensor[q_point]);
+           
+           //          if (q_point == 0)
+           //          std::cout<< stress_strain_tensor_linearized <<std::endl;
+           //          std::cout<< stress_strain_tensor <<std::endl;
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             {
+               stress_tensor = stress_strain_tensor_linearized * plast_lin_hard->get_strain(fe_values, i, q_point);
+               
+               for (unsigned int j=0; j<dofs_per_cell; ++j)
+                 {
+                   cell_matrix(i,j) += (stress_tensor *
+                                        plast_lin_hard->get_strain(fe_values, j, q_point) *
+                                        fe_values.JxW (q_point));
+                 }
+               
+               // the linearized part a(v^i;v^i,v) of the rhs
+               cell_rhs(i) += (stress_tensor *
+                               strain_tensor[q_point] *
+                               fe_values.JxW (q_point));
+               
+               // the residual part a(v^i;v) of the rhs
+               cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor *
+                               plast_lin_hard->get_strain(fe_values, i, q_point) *
+                               fe_values.JxW (q_point));
+               
+               // the residual part F(v) of the rhs
+               Tensor<1,dim> rhs_values;
+               rhs_values = 0;//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
+               //                       <<std::endl;
+             }
+         }
+       
+       for (unsigned int face=0; face<GeometryInfo<dim>::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<n_face_q_points; ++q_point)
+                 {
+                   Tensor<1,dim> 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; i<dofs_per_cell; ++i)
+                     cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
+                                     rhs_values *
+                                     fe_values_face.JxW (q_point));
+                 }
+             }
+         }
+       
+       cell->get_dof_indices (local_dof_indices);
+       constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                               local_dof_indices,
+                                               system_matrix_newton, system_rhs_newton, true);
+      };
+  
+  system_matrix_newton.compress ();
+  system_rhs_newton.compress ();
+
+  // for (unsigned int i=0; i<solution.size (); ++i)
+  // {
+  //   for (unsigned int j=0; j<solution.size (); ++j)
+  //     if (system_matrix_newton (j,i))
+  // if (constraints.is_inhomogeneously_constrained (i))
+  //     {
+  //   std::cout<< ", " <<  system_matrix_newton (j,i);
+  //   std::cout<< i << ". " << constraints.get_inhomogeneity (i)
+  //            << ", " << system_rhs_newton (i)
+  //            <<std::endl;
+  // }
+  // }
+}
 
-    std::cout << " "
-              << solver_control.last_step()
-              << " block GMRES iterations ";
-  }
+template <int dim>
+void Step4<dim>::residual_nl_system (TrilinosWrappers::MPI::Vector &u,
+                                    Vector<double>                &sigma_eff_vector)
+{
+  QGauss<dim>  quadrature_formula(2);
+  QGauss<dim-1> face_quadrature_formula(2);
+
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          UpdateFlags(update_values    |
+                                      update_gradients |
+                                      update_q_points  |
+                                      update_JxW_values));
+
+  FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula, 
+                                   update_values   | update_quadrature_points |
+                                   update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size ();
+  const unsigned int   n_face_q_points = face_quadrature_formula.size();
+
+  const EquationData::RightHandSide<dim> right_hand_side;
+  std::vector<Vector<double> > right_hand_side_values (n_q_points,
+                                                       Vector<double>(dim));  
+  std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
+                                                           Vector<double>(dim));
+  
+  Vector<double>       cell_rhs (dofs_per_cell);
+  Vector<double>       cell_sigma_eff (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  const FEValuesExtractors::Vector displacement (0);
+
+  typename DoFHandler<dim>::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<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
+       fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
+       
+       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+         {
+           SymmetricTensor<4,dim> stress_strain_tensor;
+           SymmetricTensor<2,dim> stress_tensor;
+           
+           plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point],
+                                                   elast_points, plast_points, 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 <<std::endl;*/
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             {
+               cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) *
+                               plast_lin_hard->get_strain(fe_values, i, q_point) *
+                               fe_values.JxW (q_point));
+               
+               /*          Tensor<1,dim> rhs_values;
+                           rhs_values[0] = 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<GeometryInfo<dim>::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<n_face_q_points; ++q_point)
+                 {
+                   Tensor<1,dim> 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; i<dofs_per_cell; ++i)
+                     cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
+                                     rhs_values *
+                                     fe_values_face.JxW (q_point));
+                 }
+             }
+         }
+       
+       cell->get_dof_indices (local_dof_indices);
+       constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs,
+                                                                       local_dof_indices,
+                                                                       system_rhs_newton);
+       // for (unsigned int i=0; i<dofs_per_cell; ++i)
+       //   system_rhs_newton(local_dof_indices[i]) += cell_rhs(i);
+       
+       sigma_eff_vector(cell_number) /= n_q_points;
+       cell_number += 1;
+      };
 
+  system_rhs_newton.compress ();
 
-  template <int dim>
-  void
-  StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
-  {
-    std::vector<std::string> solution_names (dim, "velocity");
-    solution_names.push_back ("pressure");
-
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-      data_component_interpretation
-      (dim, DataComponentInterpretation::component_is_part_of_vector);
-    data_component_interpretation
-      .push_back (DataComponentInterpretation::component_is_scalar);
-
-    DataOut<dim> data_out;
-    data_out.attach_dof_handler (dof_handler);
-    data_out.add_data_vector (solution, solution_names,
-                              DataOut<dim>::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 <<std::endl;
+  pcout<< "Plast-Points = " << sum_plast_points <<std::endl;
+}
 
+                                 // @sect4{Step4::projection_active_set}
 
+                                // Projection and updating of the active set
+                                 // for the dofs which penetrates the obstacle.
+template <int dim>
+void Step4<dim>::projection_active_set ()
+{
+  const EquationData::Obstacle<dim>     obstacle;
+  std::vector<bool>                     vertex_touched (dof_handler.n_dofs (), false);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end();
+  
+  TrilinosWrappers::MPI::Vector         distributed_solution (system_rhs_newton);
+  distributed_solution = solution;
+  TrilinosWrappers::MPI::Vector         lambda (solution);
+  lambda = resid_vector;
+  TrilinosWrappers::MPI::Vector         diag_mass_matrix_vector_relevant (solution);
+  diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
+
+  constraints.reinit(locally_relevant_dofs);
+  active_set.clear ();
+  const double c = 100.0*e_modul;
+  for (; cell!=endc; ++cell)
+    if (cell->is_locally_owned())
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+       if (cell->face (face)->boundary_indicator () == 9)
+         for (unsigned int v=0; v<GeometryInfo<dim-1>::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<dim> 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); 
+             
+             // std::cout<< point << ", "
+             //      << "Error: " << lambda (index_z) +
+             //   diag_mass_matrix_vector (index_z)*c*(solution_index_z - gap)
+             //      << ", " << lambda (index_z)
+             //      << ", " << diag_mass_matrix_vector (index_z)
+             //      << ", " << obstacle_value
+             //      << ", " << solution_index_z
+             //      <<std::endl;
+             
+             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_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)
+                 //           <<std::endl;
+                 // locally_relevant_dofs.print(cout);
+                 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
+                          <<std::endl;
+               }
+           }
+  
+  distributed_solution.compress(Insert);
+
+  unsigned int sum_contact_constraints = Utilities::MPI::sum(active_set.n_elements (), mpi_communicator);
+  pcout << "Number of Contact-Constaints: " << sum_contact_constraints <<std::endl;
+
+  solution = distributed_solution;
+
+  constraints.close ();
+  
+  const ConstraintMatrix::MergeConflictBehavior
+    merge_conflict_behavior = ConstraintMatrix::left_object_wins;
+  constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior);
+}
 
-  template <int dim>
-  void
-  StokesProblem<dim>::refine_mesh ()
-  {
-    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-    std::vector<bool> component_mask (dim+1, false);
-    component_mask[dim] = true;
-    KellyErrorEstimator<dim>::estimate (static_cast<const DoFHandler<dim>&>(dof_handler),
-                                        QGauss<dim-1>(degree+1),
-                                        typename FunctionMap<dim>::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 <int dim>
+void Step4<dim>::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<bool> 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<dim>(),
+                                           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<dim>(),
+                                           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<dim>(),
+                                           constraints_dirichlet_hanging_nodes,
+                                           component_mask);
+  constraints_dirichlet_hanging_nodes.close ();
+}
 
+template <int dim>
+void Step4<dim>::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<end; ++i)
+    if (constraints_hanging_nodes.is_constrained (i))
+      distributed_solution(i) = 0;
+  
+  // Solving iterative
+  SolverCG<TrilinosWrappers::MPI::Vector>
+    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() <<std::endl;
+  pcout << "   " << reduction_control.last_step()
+           << " CG iterations needed to obtain convergence with an error: "
+           <<  reduction_control.last_value()
+           << std::endl;
+
+  number_iterations += reduction_control.last_step();
+
+  constraints.distribute (distributed_solution);
+
+  solution = distributed_solution;
+}
 
+template <int dim>
+void Step4<dim>::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<std::vector<bool> > constant_modes;
+  std::vector<bool>  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 <<std::endl;     
+  // projection_active_set ();
+  IndexSet                            active_set_old (active_set);
+  Vector<double>                      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<< " " <<std::endl;
+      pcout<< j << ". Iteration of active set." <<std::endl;
+      pcout<< "Update Active Set in Dim = " << dim <<std::endl;
+      projection_active_set ();
+
+      // std::ofstream fp("Constraints");
+      // constraints.print (fp);
+
+      for (unsigned int k=0; k<=0; k++)
+       {
+         pcout<< " " <<std::endl;
+         pcout<< "Assembling ... " <<std::endl;
+         start = clock();
+         system_matrix_newton = 0;
+         system_rhs_newton = 0;
+         assemble_nl_system (solution);  //compute Newton-Matrix
+         end = clock();
+         run_time[1] += (double)(end-start)/CLOCKS_PER_SEC;
+
+         number_assemble_system += 1;
+
+         start = clock();
+         solve ();
+         end = clock();
+         run_time[2] += (double)(end-start)/CLOCKS_PER_SEC;
+
+         TrilinosWrappers::MPI::Vector    distributed_solution (system_rhs_newton);
+         distributed_solution = solution;
+
+         int damped = 0;
+         tmp_vector = old_solution;
+         double a = 0;
+         for (unsigned int i=0; (i<10)&&(!damped); i++)
+           {
+             a=pow(0.5,i);
+             old_solution = tmp_vector;
+             old_solution.sadd(1-a,a, distributed_solution);
+             
+             start = clock();
+             system_rhs_newton = 0;
+             sigma_eff_vector = 0;
+             solution = old_solution;
+             residual_nl_system (solution, sigma_eff_vector);
+             res = system_rhs_newton;
+             
+             // constraints_dirichlet_hanging_nodes.condense (res);
+             const unsigned int
+               start_res     = (res.local_range().first),
+               end_res       = (res.local_range().second);
+             for (unsigned int n=start_res; n<end_res; ++n)
+               if (constraints.is_inhomogeneously_constrained (n))
+                 {
+                   // pcout<< i << ". " << constraints.get_inhomogeneity (n)
+                   //   << ". " << res (n)
+                   //   << ", start = " << start_res
+                   //   << ", end = " << end_res
+                   //   <<std::endl;
+                   res(n) = 0;
+                 }
+
+             resid = res.l2_norm ();
+             pcout<< "Resid Newton: " << resid <<std::endl;
+             
+             if (resid<resid_old)
+               {
+                 pcout<< "--------------- alpha = " << a <<std::endl;
+                 damped=1;
+               }
+             end = clock();
+             run_time[3] = (double)(end-start)/CLOCKS_PER_SEC;
+           }
+         
+         char* name = new char[30];
+         sprintf (name,"Error-FPV.dat");
+         FILE* fp2 = fopen(name,"a");
+         
+         fprintf (fp2,"%d %le %le \n",k, a, resid);
+         fclose(fp2);
+
+         pcout<< "-------" << k << ". SQP-Iteration, error -> " << resid <<std::endl;
+
+         if (resid<1e-8)
+           {
+             pcout<< "Newton-Verfahren gestoppt bei resid = " << resid <<std::endl;
+             pcout<< "Number of Assembling systems = " << number_assemble_system <<std::endl;
+             break;
+           }
+         resid_old=resid;
+       } // End of sqp-loop
+
+      // output_results (res, "Residual");
+      resid_vector = system_rhs_newton;
+
+      if (active_set == active_set_old && resid < 1e-10)
+       break;
+      active_set_old = active_set;
+    } // End of active-set-loop
+
+  start = clock();
+  pcout<< "Creating output." <<std::endl;
+  std::ostringstream filename_solution;
+  filename_solution << "solution";
+  // filename_solution << "solution_";
+  // filename_solution << k;
+  output_results (filename_solution.str ());
+  // output_results (sigma_eff_vector, "sigma_eff");
+  end = clock();
+  run_time[4] = (double)(end-start)/CLOCKS_PER_SEC;
+
+  pcout<< "Number of Solver-Iterations = " << number_iterations <<std::endl;
+
+  pcout<< "%%%%%% Rechenzeit make grid and setup = " << run_time[0] <<std::endl;
+  pcout<< "%%%%%% Rechenzeit assemble system = " << run_time[1] <<std::endl;
+  pcout<< "%%%%%% Rechenzeit solve system = " << run_time[2] <<std::endl;
+  pcout<< "%%%%%% Rechenzeit error and lambda = " << run_time[3] <<std::endl;
+  pcout<< "%%%%%% Rechenzeit output = " << run_time[4] <<std::endl;
+}
 
-  template <int dim>
-  void StokesProblem<dim>::run ()
-  {
+template <int dim>
+void Step4<dim>::output_results (const std::string& title) const
+{
+  move_mesh (solution);
+
+  TrilinosWrappers::MPI::Vector         lambda (solution);
+  lambda = resid_vector;
+
+  DataOut<dim> data_out;
+  
+  data_out.attach_dof_handler (dof_handler);
+  
+  data_out.add_data_vector (solution, "Displacement");
+  data_out.add_data_vector (lambda, "Residual");
+  data_out.add_data_vector (active_set, "ActiveSet");
+  
+  Vector<float> subdomain (triangulation.n_active_cells());
+  for (unsigned int i=0; i<subdomain.size(); ++i)
+    subdomain(i) = triangulation.locally_owned_subdomain();
+  data_out.add_data_vector (subdomain, "subdomain");
+
+  data_out.build_patches ();
+  
+  const std::string filename = (title + "-" +
+                               Utilities::int_to_string
+                               (triangulation.locally_owned_subdomain(), 4));
+  
+  std::ofstream output_vtu ((filename + ".vtu").c_str ());
+  data_out.write_vtu (output_vtu);
+
+  if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
     {
-      std::vector<unsigned int> subdivisions (dim, 1);
-      subdivisions[0] = 1;
-
-      const Point<dim> bottom_left = (dim == 2 ?
-                                      Point<dim>(0,0) :
-                                      Point<dim>(0,0,0));
-      const Point<dim> top_right   = (dim == 2 ?
-                                      Point<dim>(1,1) :
-                                      Point<dim>(1,1,1));
-
-      GridGenerator::subdivided_hyper_rectangle (triangulation,
-                                                 subdivisions,
-                                                 bottom_left,
-                                                 top_right);
+      std::vector<std::string> filenames;
+      for (unsigned int i=0;
+          i<Utilities::MPI::n_mpi_processes(mpi_communicator);
+          ++i)
+       filenames.push_back ("solution-" +
+                            Utilities::int_to_string (i, 4) +
+                            ".vtu");
+      
+      std::ofstream master_output ((filename + ".pvtu").c_str());
+      data_out.write_pvtu_record (master_output, filenames);
     }
+  
+  TrilinosWrappers::MPI::Vector  tmp (solution);
+  tmp *= -1;
+  move_mesh (tmp);
+}
 
+template <int dim>
+void Step4<dim>::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const
+{
+  pcout<< "Moving mesh." <<std::endl;
+
+  std::vector<bool> vertex_touched (triangulation.n_vertices(),
+                                   false);
+
+  for (typename DoFHandler<dim>::active_cell_iterator
+        cell = dof_handler.begin_active ();
+       cell != dof_handler.end(); ++cell)
+    if (cell->is_locally_owned())
+      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+       {
+         if (vertex_touched[cell->vertex_index(v)] == false)
+           {
+             vertex_touched[cell->vertex_index(v)] = true;
+             
+             Point<dim> vertex_displacement;
+             for (unsigned int d=0; d<dim; ++d)
+               {
+                 if (_complete_displacement(cell->vertex_dof_index(v,d)) != 0)
+                   vertex_displacement[d]
+                     = _complete_displacement(cell->vertex_dof_index(v,d));
+               }
+             
+             cell->vertex(v) += vertex_displacement;
+           }
+       }         
+}
 
-    for (typename Triangulation<dim>::active_cell_iterator
-           cell = triangulation.begin_active();
-         cell != triangulation.end(); ++cell)
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->center()[0] == 1)
-          cell->face(f)->set_all_boundary_indicators(1);
-
-
-
-    triangulation.refine_global (1);
+template <int dim>
+void Step4<dim>::output_results (TrilinosWrappers::MPI::Vector vector,
+                                const std::string& title) const
+{
+  DataOut<dim> 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 <int dim>
+void Step4<dim>::output_results (Vector<double> vector, const std::string& title) const
+{
+  DataOut<dim> 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 <int dim>
+void Step4<dim>::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 <code>main</code> 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
+                                // <code>Step4@<2@></code> (forcing
+                                // the compiler to compile the class template
+                                // with <code>dim</code> replaced by
+                                // <code>2</code>) 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
+                                // <code>laplace_problem_2d</code> 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
+                                // <code>laplace_problem_2d</code> 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 <code>deallog.depth_console(0)</code>
+                                 // 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 <<std::endl;
+  
   return 0;
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.