]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Delete directories that have been supplanted by the single file step-35.cc.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Mar 2010 07:06:33 +0000 (07:06 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Mar 2010 07:06:33 +0000 (07:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@20885 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-35/include/EqData.h [deleted file]
deal.II/examples/step-35/include/FileReader.h [deleted file]
deal.II/examples/step-35/include/MultiComponentFunction.h [deleted file]
deal.II/examples/step-35/include/NavierStokes.h [deleted file]
deal.II/examples/step-35/source/EqData.cc [deleted file]
deal.II/examples/step-35/source/EqData_square.ccc [deleted file]
deal.II/examples/step-35/source/FileReader.cc [deleted file]
deal.II/examples/step-35/source/MultiComponentFunction.cc [deleted file]
deal.II/examples/step-35/source/NavierStokes.cc [deleted file]
deal.II/examples/step-35/source/main.cc [deleted file]

diff --git a/deal.II/examples/step-35/include/EqData.h b/deal.II/examples/step-35/include/EqData.h
deleted file mode 100644 (file)
index d528e5d..0000000
+++ /dev/null
@@ -1,71 +0,0 @@
-/*
-    Header file for the classes that define the exact solution
-    and driving force.
-
-    These classes are much like any other Function class
-    so we are not going to comment on them
-
-    by Abner Salgado.
-*/
-#ifndef _EQ_DATA_H_
-#define _EQ_DATA_H_
-
-
-/*
-    We want to be able to select which component of vector
-    valued functions we are going to work on
-*/
-#include "../include/MultiComponentFunction.h"
-
-
-// this is dealii
-#include <base/point.h>
-
-
-// We need to define sines and cosines
-#include <cmath>
-
-
-
-// This is ugly, we have to remove it
-const double PI = std::acos( -1. );
-
-
-
-// The Velocity function
-template<int dim> class Velocity: public Multi_Component_Function<dim> {
-  public:
-    Velocity(const double initial_time =0.0);
-    virtual double value(const Point<dim> &p, const unsigned int component = 0) const;
-    virtual Tensor<1,dim> gradient(const Point<dim> &p, const unsigned int component=0) const;
-    virtual void value_list( const std::vector< Point<dim> > &points, std::vector<double> &values,
-                               const unsigned int component = 0 ) const;
-    virtual void gradient_list( const std::vector< Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
-                                  const unsigned int component = 0 ) const;
-};
-
-
-
-// The Pressure function
-template<int dim> class Pressure: public Function<dim>{
-  public:
-    Pressure(const double initial_time = 0.0);
-    virtual double value(const Point<dim> &p, const unsigned int component = 0) const;
-    virtual Tensor<1,dim> gradient(const Point<dim> &p, const unsigned int component=0) const;
-    virtual void value_list( const std::vector< Point<dim> > &points, std::vector<double> &values,
-                              const unsigned int component = 0 ) const;
-    virtual void gradient_list( const std::vector< Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
-                                  const unsigned int component = 0 ) const;
-};
-
-
-
-// The Force function
-template<int dim> class Force: public Multi_Component_Function<dim>{
-  public:
-    Force( const double initial_time =0.0 );
-    virtual double value( const Point<dim> &p, const unsigned int component = 0 ) const;
-    virtual void value_list( const std::vector< Point<dim> > &points, std::vector<double> &values, const unsigned int component = 0 ) const;
-};
-
-#endif
diff --git a/deal.II/examples/step-35/include/FileReader.h b/deal.II/examples/step-35/include/FileReader.h
deleted file mode 100644 (file)
index 4047aab..0000000
+++ /dev/null
@@ -1,70 +0,0 @@
-// This class is supposed to aid in the reading of parameter files
-#ifndef _FILE_READER_H_
-#define _FILE_READER_H_
-
-
-
-#include <base/parameter_handler.h>
-
-
-
-#include <fstream>
-#include <iostream>
-
-
-
-using namespace dealii;
-
-enum Method_Formulation{
-  METHOD_STANDARD,
-  METHOD_ROTATIONAL
-};
-
-class Data_Storage{
-  public:
-    Data_Storage();
-    ~Data_Storage();
-    void read_data( char *filename );
-    void print_usage();
-    void print_status() const;
-    // The data itself
-    //// The type of method we want to use
-    Method_Formulation form;
-    //// physical data
-    double initial_time,
-           final_time,
-           Reynolds;
-    //// Time stepping data
-    double initial_dt,
-           final_dt,
-           dt_decrement;
-    //// Space discretization data
-    unsigned int n_of_global_refines,
-                 pressure_degree;
-    //// Data to solve the velocity
-    unsigned int vel_max_iterations,
-                 vel_Krylov_size,
-                 vel_off_diagonals,
-                 vel_update_prec;
-    double vel_eps,
-           vel_diag_strength;
-    //// Data to solve the projection
-    unsigned int proj_max_iterations,
-                 proj_off_diagonals;
-    double proj_eps,
-           proj_diag_strength;
-    //// Data to do the pressure update step
-    unsigned int pres_max_iterations,
-                 pres_off_diagonals;
-    double pres_eps,
-           pres_diag_strength;
-    //// Verbosity
-    bool verbose;
-    //// Frequency of the outputted data
-    unsigned int output;
-
-  protected:
-    ParameterHandler prm;
-};
-
-#endif
diff --git a/deal.II/examples/step-35/include/MultiComponentFunction.h b/deal.II/examples/step-35/include/MultiComponentFunction.h
deleted file mode 100644 (file)
index 22ba293..0000000
+++ /dev/null
@@ -1,51 +0,0 @@
-/*
-    Header file for Multi Component Functions.
-    The whole purpose of this class to exist is so we do not have to write the function
-        void set_component( const unsigned int d )
-    twice
-
-    by Abner Salgado.
-*/
-#ifndef _MULTI_COMPONENT_FUNCTION_H_
-#define _MULTI_COMPONENT_FUNCTION_H_
-
-
-/*
-  This is basically a workaround for vector valued functions
-  when we just want one of its components
-  so they are derived from the Function class
-*/
-#include <base/function.h>
-
-
-using namespace dealii;
-
-
-/*
-    The only funcitonality this class has is that it provides a common wrapper
-    for vector valued functions for the case when you only want to ask them for
-    one of their components
-*/
-template<int dim> class Multi_Component_Function: public Function<dim> {
-  public:
-    /*
-      Constructor. It does not do anything interesting but set the initial
-      time for the function to evaluate its values
-    */
-    Multi_Component_Function( const double initial_time = 0. );
-    /*
-      This is the whole reason this class exists.
-      It wouldn't seem logical to have this written twice. One for the
-      Velocity function and one for the Force function.
-      So we are instead going to derive these from this class
-      and so adding this functionality of selecting which component you
-      want to give the value of
-    */
-    void set_component(const unsigned int d );
-  protected:
-    // The current component we are working on
-    unsigned int component;
-};
-
-
-#endif
diff --git a/deal.II/examples/step-35/include/NavierStokes.h b/deal.II/examples/step-35/include/NavierStokes.h
deleted file mode 100644 (file)
index 13c8bfc..0000000
+++ /dev/null
@@ -1,278 +0,0 @@
-/*
-    Definition of the Projection Method for constant density
-    Navier Stokes
-    by Abner Salgado.
-*/
-#ifndef _NAVIER_STOKES_H_
-#define _NAVIER_STOKES_H_
-
-
-
-
-/*
-    First we include the local files, where we have the definitions of
-    1.- The exact solution
-    2.- The class that reads runtime parameters from a file
-*/
-#include "EqData.h"
-#include "FileReader.h"
-
-
-
-/*
-    These includes are to get the quadrature
-    and to handle the convergence table
-*/
-#include <base/quadrature_lib.h>
-#include <base/convergence_table.h>
-
-
-
-/*
-    The includes needed to manage multithreading
-*/
-#include <base/multithread_info.h>
-#include <base/thread_management.h>
-
-
-
-/*
-    Grid management includes.
-    Triangulation handler
-    Grid generator & refinement, accessor and iterators
-*/
-#include <grid/tria.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_refinement.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary_lib.h>
-#include <grid/grid_in.h>
-
-
-
-/*
-    Degrees of freedom management includes
-    Handling, accessing, computing useful stuff from and for them and renumbering.
-*/
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_tools.h>
-#include <dofs/dof_renumbering.h>
-#include <dofs/dof_constraints.h>
-
-
-
-/*
-    Finite element management includes
-    Definition of finite elements and extraction of info from them
-    and from finite element functions.
-*/
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_tools.h>
-
-
-
-/*
-    Linear algebra includes
-    Sparse matrices, vectors, preconditioners and solvers.
-*/
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/solver_gmres.h>
-#include <lac/sparse_ilu.h>
-#include <lac/sparse_direct.h>
-
-
-
-/*
-    Numerical algorithms includes
-    Assembly of right hand sides, computation of errors, output of the solution, etc
-*/
-#include <numerics/matrices.h>
-#include <numerics/vectors.h>
-#include <numerics/data_out.h>
-
-
-/*
-    This class is the core of all the program
-    It is the implementation of the projection method for the Navier-Stokes equations.
-*/
-template<int dim> class Navier_Stokes_Projection{
-  public:
-    /*
-      Constructor. It takes as a parameter a reference to a Data_Storage class, which is
-      basically a container of all the data this class needs to be defined (and some other stuff)
-      Here we basically just copy the needed data and/or assign the values that are needed.
-    */
-    Navier_Stokes_Projection( const Data_Storage &data );
-    /*
-      Destructor. Completely trivial.
-    */
-    ~Navier_Stokes_Projection();
-    /*
-      This function creates the triangulation and refines it the required number of times
-      It also initializes all the quantities that are mesh dependent but not time stepping dependent
-      i.e. the size of matrices and vectors.
-    */
-    void Create_Triangulation( const unsigned int n_of_refines );
-    /*
-      Having created a mesh and once the time step is set, the next step is to initialize the method
-      i.e. compute the matrices that are never going to change, load the initial data, etc.
-      This is what this function does.
-    */
-    void Initialize();
-    /*
-      This is the time marching function, which starting at t_0 advances in time using the projection method
-      with time step dt until T.
-      The boolean parameter that it takes is to enable information about what the function is doing at the present
-      moment, i.e. diffusion, projection substep; updating preconditioners etc.
-      This is useful mostly for debuggin purposes and so it is by default set to false
-    */
-    void run( const bool verbose = false, const unsigned int n_of_plots = 10 );
-    /*
-      Having reached the final time T, we want to measure the error that we have made.
-      This method is responsible for that. Saves the results in a ConvergenceTable object
-      which later we can print or compute things with it
-    */
-    void Post_Process();
-    /*
-      The whole reason for this class, and all this code after all, is to make convergence tests
-      with respect to the time discretization. For this reason we need to be able to vary the time step
-      we are going to work with. This method sets up the time step that is to be used in the given run
-    */
-    void set_dt( const double ddt );
-
-  protected:
-    // The type of method
-    Method_Formulation type;
-
-    // Discretization data
-    //// The polynomial degree
-    unsigned int deg;
-    //// The time step
-    double dt;
-
-    // Physical data:
-    //// initial time, final time, Reynolds number
-    double t_0, T, Re;
-    //// The external driving force
-    Force<dim> rhs;
-    //// The exact velocity (to apply boundary values)
-    Velocity<dim> vel_exact;
-    //// The boundary conditions
-    std::map<unsigned int, double> boundary_values;
-
-    // Finite Element Spaces data
-    //// The mesh
-    Triangulation<dim> triangulation;
-    //// The DoF handlers
-    DoFHandler<dim> dof_handler_velocity, dof_handler_pressure;
-    //// The polynomial spaces
-    FE_Q<dim> fe_velocity, fe_pressure;
-    //// The quadrature formulae
-    QGauss<dim> quadrature_pressure, quadrature_velocity;
-
-    /*
-      Linear Algebra Data
-        The sparsity patterns where the matrices will live
-    */
-    SparsityPattern spar_pattern_velocity, spar_pattern_pressure, spar_pattern_pres_vel;
-    /*
-      The actual matrices. The projection matrix never changes.
-      For the velocity a part of the matrix never changes (neither in time nor for each component),
-      namely the part related with the time derivative and the diffusion
-      but, the advection part changes and so we need dim+1 matrices. One to store
-      the constant part, and dim for each matrix at each iteration
-    */
-    SparseMatrix<double> vel_Laplace_plus_Mass, vel_it_matrix[dim], vel_Mass, vel_Laplace,
-                         pres_Laplace, pres_Mass, pres_Diff[dim];
-    /*
-      We need to regularize the Laplace operator for the pressure
-      for that we use a constraint
-    */
-    ConstraintMatrix pres_regularization;
-    //// The solutions at times n and (n-1)
-    Vector<double> pres_n, pres_n_minus_1, phi_n, phi_n_minus_1, u_n[dim], u_n_minus_1[dim],
-    //// The time extrapolation of the velocity, used to make the semi-implicit time discretization
-    //// of the advection term
-                  u_star[dim],
-    //// Right hand side for the component of the momentum equation
-                  force[dim],
-    //// Temporary arrays
-    v_tmp, pres_tmp;
-    //// The preconditioners
-//     SparseILU<double> prec_velocity[dim];
-    SparseDirectUMFPACK prec_velocity[dim], prec_mass, prec_pressure;
-
-    // The convergence table
-    ConvergenceTable convergence_table;
-
-    // Exception we will throw if the time step is invalid
-    DeclException2( ExcInvalidTimeStep, double, double, <<" The time step "<<arg1<<" is out of range."<<std::endl
-                                                                <<" The permitted range is (0,"<<arg2<<"]");
-
-    /*
-      A diffusion step. The boolean parameter that it takes indicates if we want the preconditioners
-      for each component of the velocity to be updated. Recall that the operator that these matrices
-      represent is of the form
-          1.5/dt Mass + 1/Re Laplace_h + ( u^*.D )
-      so the dominating part is
-          1.5/dt Mass + 1/Re Laplace_h
-      and we can recycle the preconditioner for several steps to save computational time
-    */
-    inline void diffusion_step( const bool reinit_prec );
-    /*
-        A projection step. The parameter that it takes indicates if we want the preconditioner to be
-        updated. Since this matrix never changes, ideally this should be done just once.
-    */
-    inline void projection_step( const bool reinit_prec );
-    /*
-      A pressure update step.
-      Capable of recognizing which formulation of the method we are using.
-      If it is the standard formulation, then the pressure update is done by
-          p^{n+1} = p^{n} + phi^{n+1}
-      which involves just adding numerical values.
-      On the other hand, for the rotational formulation of the method the pressure update
-      is done by the formula
-          p^{n+1} = p^{n} + phi^{n+1} +1/Re div u^{n+1}
-      which in any case involves the solution of a mass matrix problem.
-      NOTE: Here we assume that the mesh is uniformly refined, and thus we do not use
-      a preconditioner for this mass matrix problem
-    */
-    inline void update_pressure( const bool reinit_prec );
-  private:
-
-    // A bunch of data for the solvers and preconditioners
-    // data solve velocity
-    unsigned int vel_max_its, vel_Krylov_size, vel_off_diagonals, vel_update_prec;
-    double vel_eps, vel_diag_strength;
-    // data solve projection
-    unsigned int proj_max_its, proj_off_diagonals;
-    double proj_eps, proj_diag_strength;
-    // data solve pressure update
-    unsigned int pres_max_its, pres_off_diagonals;
-    double pres_eps, pres_diag_strength;
-
-    // Extremely simple function that computes the velocity extrapolant. Added for readability
-    inline void interpolate_velocity();
-
-    //////// initialization of constant wrt dt matrices
-    inline void init_velocity_matrices();
-    inline void init_pressure_matrices();
-    inline void init_gradient_operator();
-
-    // Assembly of the advection term
-    inline void assemble_advection_term( const unsigned int d );
-
-    // Solution of each component of the velocity on the diffusion step
-    inline void threaded_solve( const unsigned int d );
-    inline void plot_solution( const unsigned int step );
-};
-
-
-
-#endif
diff --git a/deal.II/examples/step-35/source/EqData.cc b/deal.II/examples/step-35/source/EqData.cc
deleted file mode 100644 (file)
index 8a8e198..0000000
+++ /dev/null
@@ -1,170 +0,0 @@
-/*
-    Implementation of the classes that represent the exact solution
-    and external force FOR THE PROBLEM ON THE DISK
-
-
-    There's not much to say here. All these classes just represent
-    mathematical formulae
-
-    by Abner Salgado.
-*/
-
-
-
-#include "../include/EqData.h"
-
-
-
-// Force function methods
-template<int dim> Force<dim>::Force( const double initial_time ): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void  Force<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
-                                                  const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    values[i] = Force<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Force<dim>::value(const Point<dim> &p, const unsigned int) const{
-  double t = FunctionTime::get_time(),
-         cosx = std::cos( p(0) ),
-         sinx = std::sin( p(0) ),
-         cosy = std::cos( p(1) ),
-         siny = std::sin( p(1) ),
-         cost = std::cos(t),
-         sint = std::sin(t),
-         return_value = 0.;
-  switch( Multi_Component_Function<dim>::component ){
-    case 0:
-    // y*sin(t)-x*cos(t)^2+cos(x)*sin(y)*sin(t)
-      return_value = p(1)*sint - p(0)*cost*cost + cosx*siny*sint;
-      break;
-    case 1:
-    // -x*sin(t)-y*cos(t)^2+sin(x)*cos(y)*sin(t)
-      return_value = -p(0)*sint - p(1)*cost*cost + sinx*cosy*sint ;
-
-      break;
-    default:
-      Assert( false, ExcNotImplemented() );
-  };
-  return return_value;
-}
-
-
-
-// Velocity function methods
-template<int dim> Velocity<dim>::Velocity(const double initial_time): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void  Velocity<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
-                                                  const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    values[i] = Velocity<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Velocity<dim>::value(const Point<dim> &p, const unsigned int) const{
-  double return_value = std::cos( Function<dim>::get_time() );
-  switch( Multi_Component_Function<dim>::component ){
-    case 0:
-     // -y*cos(t)
-       return_value *= -p(1);
-      break;
-    case 1:
-      // x*cos(t)
-      return_value *= p(0);
-      break;
-    default:
-      Assert( false, ExcNotImplemented() );
-  };
-  return return_value;
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Velocity<dim>::gradient(const Point<dim> &p, const unsigned int) const{
-  Tensor<1,dim> return_value;
-  switch( Multi_Component_Function<dim>::component ){
-    // [0, -cos(t)]
-    case 0:
-     return_value[0] = 0.;
-     return_value[1] = -std::cos( Function<dim>::get_time() );
-      break;
-    case 1:
-    // [cos(t), 0]
-     return_value[0] = std::cos( Function<dim>::get_time() );
-     return_value[1] = 0.;
-      break;
-    default:
-      Assert( false, ExcNotImplemented() );
-  };
-  return return_value;
-}
-
-
-
-template<int dim> void Velocity<dim>::gradient_list( const std::vector<Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
-                                                        const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( gradients.size() == n_points, ExcDimensionMismatch( gradients.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    gradients[i] = Velocity<dim>::gradient( points[i] );
-}
-
-
-
-// Pressure function methods
-template<int dim> Pressure<dim>::Pressure(const double initial_time): Function<dim>(1,initial_time){}
-
-
-
-template<int dim> inline double Pressure<dim>::value(const Point<dim> &p, const unsigned int) const{
-  // sin(x-y+t)
-  return std::sin( p(0) )*std::sin( p(1) )*std::sin( Function<dim>::get_time() );
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Pressure<dim>::gradient(const Point<dim> &p, const unsigned int) const{
-  // [cos(x)*sin(y)*sin(t), sin(x)*cos(y)*sin(t)]
-  return Point<dim>( std::cos( p(0) )*std::sin( p(1) )*std::sin( Function<dim>::get_time() ),
-                       std::sin( p(0) )*std::cos( p(1) )*std::sin( Function<dim>::get_time() ) );
-}
-
-
-
-template<int dim> void Pressure<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
-                                                    const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    values[i] = Pressure<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline void Pressure<dim>::gradient_list( const std::vector<Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
-                                                                const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( gradients.size() == n_points, ExcDimensionMismatch( gradients.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    gradients[i] = Pressure<dim>::gradient( points[i] );
-}
-
-
-
-// explicit template instantiation
-template class Force<deal_II_dimension>;
-template class Velocity<deal_II_dimension>;
-template class Pressure<deal_II_dimension>;
diff --git a/deal.II/examples/step-35/source/EqData_square.ccc b/deal.II/examples/step-35/source/EqData_square.ccc
deleted file mode 100644 (file)
index 2f52b5a..0000000
+++ /dev/null
@@ -1,211 +0,0 @@
-/*
-    Implementation of the classes that represent the exact solution
-    and external force.
-
-    There's not much to say here. All these classes just represent
-    mathematical formulae
-
-    by Abner Salgado.
-*/
-
-
-
-#include "../include/EqData.h"
-
-
-
-// Force function methods
-template<int dim> Force<dim>::Force( const double initial_time ): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void  Force<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
-                                                  const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    values[i] = Force<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Force<dim>::value(const Point<dim> &p, const unsigned int) const{
-  double t = FunctionTime::get_time(),
-         sin_2pi_x  = std::sin( 2.*PI*p(0) ),
-         sin_2pi_y  = std::sin( 2.*PI*p(1) ),
-         cos_2pi_x  = std::cos( 2.*PI*p(0) ),
-         cos_2pi_y  = std::cos( 2.*PI*p(1) ),
-         sin_pi_x   = std::sin( PI*p(0) ),
-         sin_pi_y   = std::sin( PI*p(1) ),
-         cos_pi_x   = std::cos( PI*p(0) ),
-         cos_pi_y   = std::cos( PI*p(1) ),
-         cos_t      = std::cos( t ),
-         sin_t      = std::sin( t ),
-         return_value = 0.;
-  switch( Multi_Component_Function<dim>::component ){
-    case 0:
-    // Pi*sin(2*Pi*y)*sin(Pi*x)^2*cos(t)
-    // -2*Pi^3*sin(2*Pi*y)*cos(Pi*x)^2*sin(t)
-    // +6*Pi^3*sin(2*Pi*y)*sin(Pi*x)^2*sin(t)
-    // +2*Pi^3*sin(2*Pi*y)^2*sin(Pi*x)^3*sin(t)^2*cos(Pi*x)
-    // -2*Pi^3*sin(2*Pi*x)*sin(Pi*y)^2*sin(t)^2*cos(2*Pi*y)*sin(Pi*x)^2
-    // -Pi*sin(Pi*x)*sin(Pi*y)*sin(t)
-      return_value = PI*sin_2pi_y*sin_pi_x*sin_pi_x*cos_t
-                   - 2.*PI*PI*PI*sin_2pi_y*cos_pi_x*cos_pi_x*sin_t
-                   + 6.*PI*PI*PI*sin_2pi_y*sin_pi_x*sin_pi_x*sin_t
-                   + 2.*PI*PI*PI*sin_2pi_y*sin_2pi_y*sin_pi_x*sin_pi_x*sin_pi_x*sin_t*sin_t*cos_pi_x
-                   - 2.*PI*PI*PI*sin_2pi_x*sin_pi_y*sin_pi_y*sin_t*sin_t*cos_2pi_y*sin_pi_x*sin_pi_x
-                   - PI*sin_pi_x*sin_pi_y*sin_t;
-
-      break;
-    case 1:
-    // -Pi*sin(2*Pi*x)*sin(Pi*y)^2*cos(t)
-    // -6*Pi^3*sin(2*Pi*x)*sin(Pi*y)^2*sin(t)
-    // +2*Pi^3*sin(2*Pi*x)*cos(Pi*y)^2*sin(t)
-    // +Pi*cos(Pi*x)*cos(Pi*y)*sin(t)
-    // -2*Pi^3*sin(2*Pi*y)*sin(Pi*x)^2*sin(t)^2*cos(2*Pi*x)*sin(Pi*y)^2
-    // +2*Pi^3*sin(2*Pi*x)^2*sin(Pi*y)^3*sin(t)^2*cos(Pi*y)
-      return_value = -PI*sin_2pi_x*sin_pi_y*sin_pi_y*cos_t
-                   - 6.*PI*PI*PI*sin_2pi_x*sin_pi_y*sin_pi_y*sin_t
-                   + 2.*PI*PI*PI*sin_2pi_x*cos_pi_y*cos_pi_y*sin_t
-                   + PI*cos_pi_x*cos_pi_y*sin_t
-                   - 2.*PI*PI*PI*sin_2pi_y*sin_pi_x*sin_pi_x*sin_t*sin_t*cos_2pi_x*sin_pi_y*sin_pi_y
-                   + 2.*PI*PI*PI*sin_2pi_x*sin_2pi_x*sin_pi_y*sin_pi_y*sin_pi_y*sin_t*sin_t*cos_pi_y;
-
-      break;
-    default:
-      Assert( false, ExcNotImplemented() );
-  };
-  return return_value;
-}
-
-
-
-// Velocity function methods
-template<int dim> Velocity<dim>::Velocity(const double initial_time): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void  Velocity<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
-                                                  const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    values[i] = Velocity<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Velocity<dim>::value(const Point<dim> &p, const unsigned int) const{
-  double sin_2pi_x  = std::sin( 2.*PI*p(0) ),
-         sin_2pi_y  = std::sin( 2.*PI*p(1) ),
-         sin_pi_x   = std::sin( PI*p(0) ),
-         sin_pi_y   = std::sin( PI*p(1) ),
-         sin_t      = std::sin( Function<dim>::get_time() ),
-         return_value = 0.;
-  switch( Multi_Component_Function<dim>::component ){
-    case 0:
-     // pi*sin(2*pi*y)*sin(pi*x)^2*sin(t)
-       return_value = PI*sin_2pi_y*sin_pi_x*sin_pi_x*sin_t;
-      break;
-    case 1:
-      // -pi*sin(2*pi*x)*sin(pi*y)^2*sin(t)
-      return_value = -PI*sin_2pi_x*sin_pi_y*sin_pi_y*sin_t;
-      break;
-    default:
-      Assert( false, ExcNotImplemented() );
-  };
-  return return_value;
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Velocity<dim>::gradient(const Point<dim> &p, const unsigned int) const{
-  Tensor<1,dim> return_value;
-  double sin_2pi_x = std::sin( 2.*PI*p(0) ),
-         sin_2pi_y = std::sin( 2.*PI*p(1) ),
-         cos_2pi_x = std::cos( 2.*PI*p(0) ),
-         cos_2pi_y = std::cos( 2.*PI*p(1) ),
-         sin_pi_x  = std::sin( PI*p(0) ),
-         sin_pi_y  = std::sin( PI*p(1) ),
-         cos_pi_x  = std::cos( PI*p(0) ),
-         cos_pi_y  = std::cos( PI*p(1) ),
-         sin_t     = std::sin( Function<dim>::get_time() );
-  switch( Multi_Component_Function<dim>::component ){
-    // [2*Pi^2*sin(2*Pi*y)*sin(Pi*x)*sin(t)*cos(Pi*x),
-    //     2*Pi^2*cos(2*Pi*y)*sin(Pi*x)^2*sin(t)]
-    case 0:
-     return_value[0] = 2.*PI*PI*sin_2pi_y*sin_pi_x*sin_t*cos_pi_x;
-     return_value[1] = 2.*PI*PI*cos_2pi_y*sin_pi_x*sin_pi_x*sin_t;
-      break;
-    case 1:
-    // [-2*Pi^2*cos(2*Pi*x)*sin(Pi*y)^2*sin(t),
-    // -2*Pi^2*sin(2*Pi*x)*sin(Pi*y)*sin(t)*cos(Pi*y)]
-     return_value[0] = -2.*PI*PI*cos_2pi_x*sin_pi_y*sin_pi_y*sin_t;
-     return_value[1] = -2.*PI*PI*sin_2pi_x*sin_pi_y*sin_t*cos_pi_y;
-      break;
-    default:
-      Assert( false, ExcNotImplemented() );
-  };
-  return return_value;
-}
-
-
-
-template<int dim> void Velocity<dim>::gradient_list( const std::vector<Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
-                                                        const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( gradients.size() == n_points, ExcDimensionMismatch( gradients.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    gradients[i] = Velocity<dim>::gradient( points[i] );
-}
-
-
-
-// Pressure function methods
-template<int dim> Pressure<dim>::Pressure(const double initial_time): Function<dim>(1,initial_time){}
-
-
-
-template<int dim> inline double Pressure<dim>::value(const Point<dim> &p, const unsigned int) const{
-  // cos(pi*x)*sin(pi*y)*sin(t)
-  return std::cos( PI*p(0) )*std::sin( PI*p(1) )*std::sin( Function<dim>::get_time() );
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Pressure<dim>::gradient(const Point<dim> &p, const unsigned int) const{
-  // || -sin(Pi*x)*Pi*sin(Pi*y)*sin(t) ||
-  // || cos(Pi*x)*cos(Pi*y)*Pi*sin(t)  ||
-  return ( PI*std::sin( Function<dim>::get_time() ) )* Point<dim>( -std::sin( PI*p(0) )*std::sin( PI*p(1) ),
-                                                                    std::cos( PI*p(0) )*std::cos( PI*p(1) ) ) ;
-}
-
-
-
-template<int dim> void Pressure<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
-                                                    const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    values[i] = Pressure<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline void Pressure<dim>::gradient_list( const std::vector<Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
-                                                                const unsigned int ) const{
-  const unsigned int n_points = points.size();
-  Assert( gradients.size() == n_points, ExcDimensionMismatch( gradients.size(), n_points ) );
-  for (unsigned int i=0; i<n_points; ++i)
-    gradients[i] = Pressure<dim>::gradient( points[i] );
-}
-
-
-
-// explicit template instantiation
-template class Force<deal_II_dimension>;
-template class Velocity<deal_II_dimension>;
-template class Pressure<deal_II_dimension>;
diff --git a/deal.II/examples/step-35/source/FileReader.cc b/deal.II/examples/step-35/source/FileReader.cc
deleted file mode 100644 (file)
index ce58a8c..0000000
+++ /dev/null
@@ -1,187 +0,0 @@
-#include "../include/FileReader.h"
-
-
-
-// Constructor. Here we set up the
-// Data Format we are going to use
-Data_Storage::Data_Storage(){
-  prm.declare_entry( "Method_Form", "rotational", Patterns::Selection( "rotational|standard" ),
-                    " Used to select the type of method that we are going to use. " );
-
-  // Physical data of the problem
-  prm.enter_subsection( "Physical data" );
-    prm.declare_entry( "initial_time", "0.", Patterns::Double( 0. ), " The initial time of the simulation. " );
-    prm.declare_entry( "final_time", "1.", Patterns::Double( 0. ), " The final time of the simulation. " );
-    prm.declare_entry( "Reynolds", "1.", Patterns::Double( 0. ), " The Reynolds number. " );
-  prm.leave_subsection();
-
-  // Time stepping data of the problem
-  prm.enter_subsection( "Time step data" );
-    prm.declare_entry( "initial_dt", "0.1", Patterns::Double( 0. ), " The initial time step size. " );
-    prm.declare_entry( "final_dt", "5e-4",  Patterns::Double( 0. ), " The final time step size. " );
-    prm.declare_entry( "dt_decrement", "2.", Patterns::Double( 1.5 ), " The factor by which the time step will be divided. " );
-  prm.leave_subsection();
-
-  // Space discretization data
-  prm.enter_subsection( "Space discretization" );
-    prm.declare_entry( "n_of_refines", "5", Patterns::Integer( 1, 15), " The number of global refines we do on the mesh. " );
-    prm.declare_entry( "pressure_fe_degree", "1", Patterns::Integer( 1, 5 ), " The polynomial degree for the pressure space. " );
-  prm.leave_subsection();
-
-  // Velocity solution data
-  prm.enter_subsection( "Data solve velocity" );
-    prm.declare_entry( "max_iterations", "1000", Patterns::Integer( 1, 1000 ), " The maximal number of iterations GMRES must make. " );
-    prm.declare_entry( "eps", "1e-12", Patterns::Double( 0. ), " The stopping criterion. " );
-    prm.declare_entry( "Krylov_size", "30", Patterns::Integer(1), " The size of the Krylov subspace to be used. " );
-    prm.declare_entry( "off_diagonals", "60", Patterns::Integer(0), " The number of off-diagonal elements ILU must compute. " );
-    prm.declare_entry( "diag_strength", "0.01", Patterns::Double( 0. ), " Diagonal strengthening coefficient. " );
-    prm.declare_entry( "update_prec", "15", Patterns::Integer(1), " This number indicates how often we need to update the preconditioner" );
-  prm.leave_subsection();
-
-  // Projection step data
-  prm.enter_subsection( "Data solve projection" );
-    prm.declare_entry( "max_iterations", "1000", Patterns::Integer( 1, 1000 ), " The maximal number of iterations CG must make. " );
-    prm.declare_entry( "eps", "1e-12", Patterns::Double( 0. ), " The stopping criterion. " );
-    prm.declare_entry( "off_diagonals", "100", Patterns::Integer(1), " The number of off-diagonal elements ILU must compute" );
-    prm.declare_entry( "diag_strength", "0.1", Patterns::Double( 0. ), " Diagonal strengthening coefficient. " );
-  prm.leave_subsection();
-
-  // Pressure update data
-  prm.enter_subsection( "Data solve pressure update" );
-    prm.declare_entry( "max_iterations", "1000", Patterns::Integer( 1, 1000 ), " The maximal number of iterations CG must make. " );
-    prm.declare_entry( "eps", "1e-12", Patterns::Double( 0. ), " The stopping criterion. " );
-    prm.declare_entry( "off_diagonals", "10", Patterns::Integer(0), " The number of off-diagonal elements that ILU must compute" );
-    prm.declare_entry( "diag_strength", "0.", Patterns::Double(0), " Diagonal strengthening coefficient" );
-  prm.leave_subsection();
-
-  // Verbosity of output
-  prm.declare_entry( "verbose", "true", Patterns::Bool(), " This indicates whether the output of the solution process should be verbose. " );
-
-  // How often we want the data to be outputted
-  prm.declare_entry( "output", "10", Patterns::Integer(1), " This indicates between how many time steps we print the solution. " );
-}
-
-
-
-// Destructor. Does nothing
-Data_Storage::~Data_Storage(){
-}
-
-
-
-// Here is where all happens. We read all the data from the indicated file
-void Data_Storage::read_data( char *filename ){
-  std::ifstream file;
-  file.open( filename );
-  if( not file )
-    throw ExcFileNotOpen( filename );
-  prm.read_input( file );
-
-  std::string token = prm.get( "Method_Form" );
-  if( token == std::string("rotational") )
-    form = METHOD_ROTATIONAL;
-  else
-    form = METHOD_STANDARD;
-
-  // Physical data of the problem
-  prm.enter_subsection( "Physical data" );
-    initial_time = prm.get_double( "initial_time" );
-    final_time   = prm.get_double( "final_time" );
-    Reynolds     = prm.get_double( "Reynolds" );
-  prm.leave_subsection();
-
-  // Time stepping data of the problem
-  prm.enter_subsection( "Time step data" );
-    initial_dt   = prm.get_double( "initial_dt" );
-    final_dt     = prm.get_double( "final_dt" );
-    dt_decrement = prm.get_double( "dt_decrement" );
-  prm.leave_subsection();
-
-  // Space discretization data
-  prm.enter_subsection( "Space discretization" );
-    n_of_global_refines = prm.get_integer( "n_of_refines" );
-    pressure_degree     = prm.get_integer( "pressure_fe_degree" );
-  prm.leave_subsection();
-
-  // Velocity solution data
-  prm.enter_subsection( "Data solve velocity" );
-    vel_max_iterations = prm.get_double( "max_iterations" );
-    vel_eps            = prm.get_double( "eps" );
-    vel_Krylov_size    = prm.get_integer( "Krylov_size" );
-    vel_off_diagonals  = prm.get_integer( "off_diagonals" );
-    vel_diag_strength  = prm.get_double( "diag_strength" );
-    vel_update_prec    = prm.get_integer( "update_prec" );
-  prm.leave_subsection();
-
-  // Projection step data
-  prm.enter_subsection( "Data solve projection" );
-    proj_max_iterations = prm.get_integer( "max_iterations" );
-    proj_eps            = prm.get_double( "eps" );
-    proj_off_diagonals  = prm.get_integer( "off_diagonals" );
-    proj_diag_strength  = prm.get_double( "diag_strength" );
-  prm.leave_subsection();
-
-  // Pressure update data
-  prm.enter_subsection( "Data solve pressure update" );
-    pres_max_iterations = prm.get_integer( "max_iterations" );
-    pres_eps            = prm.get_double( "eps" );
-    pres_off_diagonals  = prm.get_integer( "off_diagonals" );
-    pres_diag_strength  = prm.get_double( "diag_strength" );
-  prm.leave_subsection();
-
-  // Verbosity
-  verbose = prm.get_bool( "verbose" );
-
-  // Output frequency
-  output = prm.get_integer( "output" );
-
-  file.close();
-}
-
-
-
-// Prints the current values of the data
-// Mostly (if not only) used for debugging purposes
-void Data_Storage::print_status() const{
-  std::cout<<"Method Form = "<<( (form == METHOD_ROTATIONAL)?"rotational":"standard" )<<std::endl
-           <<"Physical Data:"<<std::endl
-           <<"  initial_time = "<<initial_time<<std::endl
-           <<"  final_time = "<<final_time<<std::endl
-           <<"  Re = "<<Reynolds<<std::endl<<std::endl
-           <<"Time step data:"<<std::endl
-           <<"  initial_dt = "<<initial_dt<<std::endl
-           <<"  final_dt = "<<final_dt<<std::endl
-           <<"  dt_decrement = "<<dt_decrement<<std::endl<<std::endl
-           <<"Space discretization data:"<<std::endl
-           <<"  n_of_global_refines = "<<n_of_global_refines<<std::endl
-           <<"  Pressure space degree = "<<pressure_degree<<std::endl<<std::endl
-           <<"Data for the solution of the diffusion:"<<std::endl
-           <<"  max iterations = "<<vel_max_iterations<<std::endl
-           <<"  eps = "<<vel_eps<<std::endl
-           <<"  Krylov subspace size = "<<vel_Krylov_size<<std::endl
-           <<"  off diagonals = "<<vel_off_diagonals<<std::endl
-           <<"  diagonal strengthening = "<<vel_diag_strength<<std::endl
-           <<"  update prec every "<<vel_update_prec<<" steps"<<std::endl<<std::endl
-           <<"Data for the projection step:"<<std::endl
-           <<"  max iterations = "<<proj_max_iterations<<std::endl
-           <<"  eps = "<<proj_eps<<std::endl
-           <<"  off diagonals = "<<proj_off_diagonals<<std::endl
-           <<"  diagonal strengthening = "<<proj_diag_strength<<std::endl<<std::endl
-           <<"Data for the pressure update step"<<std::endl
-           <<"  max iterations = "<<pres_max_iterations<<std::endl
-           <<"  eps = "<<pres_eps<<std::endl
-           <<"  off_diagonals = "<<pres_off_diagonals<<std::endl
-           <<"  diagonal strengthening = "<<pres_diag_strength<<std::endl<<std::endl
-           <<"Verbose = "<<(verbose?"true":"false")<<std::endl<<std::endl
-           <<"Output ="<<output<<std::endl<<std::endl;
-}
-
-
-
-// This method just prints the format the
-// input data file should have
-void Data_Storage::print_usage(){
-  static const char *message = "";
-  std::cout<<message<<std::endl;
-  prm.print_parameters( std::cout, ParameterHandler::Text );
-}
diff --git a/deal.II/examples/step-35/source/MultiComponentFunction.cc b/deal.II/examples/step-35/source/MultiComponentFunction.cc
deleted file mode 100644 (file)
index 5624c9d..0000000
+++ /dev/null
@@ -1,26 +0,0 @@
-/*
-    Muti_Component_Function class implementation
-*/
-#include "../include/MultiComponentFunction.h"
-
-
-// Constructor: Just call the Function constructor and set the component to 0
-template<int dim> Multi_Component_Function<dim>::Multi_Component_Function( const double initial_time ):
-                    Function<dim>( 1, initial_time ), component(0) {
-}
-
-
-
-// Set Component Function: Check that it is in range and then set it
-template<int dim> void Multi_Component_Function<dim>::set_component(const unsigned int d ){
-  // Check if the requested component is correct
-  Assert( d<dim, ExcIndexRange( d, 0, dim ) );
-  // We have only written 2d
-  Assert( dim >= 2, ExcNotImplemented() );
-  component = d;
-}
-
-
-
-// explicit template instantiation
-template class Multi_Component_Function<deal_II_dimension>;
diff --git a/deal.II/examples/step-35/source/NavierStokes.cc b/deal.II/examples/step-35/source/NavierStokes.cc
deleted file mode 100644 (file)
index 88961bd..0000000
+++ /dev/null
@@ -1,792 +0,0 @@
-/*
-    Implementation of the Navier_Stokes_projection class
-    by Abner Salgado.
-*/
-#include "../include/NavierStokes.h"
-
-
-
-#include <sys/time.h>
-
-/// debug
-bool show;
-double Solve_time;
-///----
-
-
-// Constructor
-template<int dim> Navier_Stokes_Projection<dim>::Navier_Stokes_Projection( const Data_Storage &data):
-                    type( data.form ), deg( data.pressure_degree ), dt( data.initial_dt ), t_0( data.initial_time ),
-                    T( data.final_time ), Re( data.Reynolds ), rhs( data.initial_time ), vel_exact( data.initial_time ),
-                    dof_handler_velocity(triangulation), dof_handler_pressure(triangulation),
-                    fe_velocity(deg+1), fe_pressure(deg), quadrature_pressure(deg+1), quadrature_velocity(deg+2),
-                    vel_max_its( data.vel_max_iterations ), vel_Krylov_size( data.vel_Krylov_size ),
-                    vel_off_diagonals( data.vel_off_diagonals ),
-                    vel_update_prec( data.vel_update_prec ), vel_eps( data.vel_eps ), vel_diag_strength( data.vel_diag_strength),
-                    proj_max_its( data.proj_max_iterations ), proj_off_diagonals( data.proj_off_diagonals ),
-                    proj_eps( data.proj_eps ), proj_diag_strength( data.proj_diag_strength ),
-                    pres_max_its( data.pres_max_iterations), pres_off_diagonals( data.pres_off_diagonals ),
-                    pres_eps( data.pres_eps ), pres_diag_strength( data.pres_diag_strength )
-{
-  // After having initialized the bunch of data we do nothing
-  // NOTE TO SELF: Do I need to do this check?
-  if(deg < 1)
-  std::cout<<" WARNING: The chosen pair of finite element spaces is not stable."<<std::endl
-           <<" The obtained results will be nonsense"<<std::endl;
-  // We check that the time step is within the allowed limits
-  AssertThrow( not ( ( dt <= 0. ) or ( dt > .5*T ) ), ExcInvalidTimeStep( dt, .5*T ) );
-}
-
-
-
-// Destructor
-template<int dim> Navier_Stokes_Projection<dim>::~Navier_Stokes_Projection(){
-  dof_handler_velocity.clear();
-  dof_handler_pressure.clear();
-}
-
-
-
-// Set time step
-template<int dim> void Navier_Stokes_Projection<dim>::set_dt( const double ddt ){
-  // We just check that it is within the permitted limits
-  AssertThrow( not ( ( ddt <= 0. ) or ( ddt > .5*T ) ), ExcInvalidTimeStep( ddt, .5*T ) );
-  dt = ddt;
-}
-
-
-
-// Initialization of the velocity matrices and assembly of those that do not depend on dt
-template<int dim> void Navier_Stokes_Projection<dim>::init_velocity_matrices(){
-  //// Init the sparsity pattern for the velocity
-  spar_pattern_velocity.reinit( dof_handler_velocity.n_dofs(), dof_handler_velocity.n_dofs(),
-                                dof_handler_velocity.max_couplings_between_dofs() );
-  DoFTools::make_sparsity_pattern( dof_handler_velocity, spar_pattern_velocity );
-  spar_pattern_velocity.compress();
-
-  //// Init the matrices for the velocity
-  vel_Laplace_plus_Mass.reinit( spar_pattern_velocity );
-  for( unsigned int d=0; d<dim; ++d )
-    vel_it_matrix[d].reinit( spar_pattern_velocity );
-  vel_Mass.reinit( spar_pattern_velocity );
-  vel_Laplace.reinit( spar_pattern_velocity );
-
-  // We assemble the mass matrix
-  MatrixCreator::create_mass_matrix( dof_handler_velocity, quadrature_velocity, vel_Mass );
-
-  // We assemble the Laplace matrix
-  MatrixCreator::create_laplace_matrix( dof_handler_velocity, quadrature_velocity, vel_Laplace );
-}
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::init_pressure_matrices(){
-  //// Init the sparsity pattern for the pressure
-  spar_pattern_pressure.reinit( dof_handler_pressure.n_dofs(), dof_handler_pressure.n_dofs(),
-                                dof_handler_pressure.max_couplings_between_dofs() );
-  DoFTools::make_sparsity_pattern( dof_handler_pressure, spar_pattern_pressure );
-
-  // Before we close the sparsity pattern, we need to
-  // init the constraints for the Laplace operator on the pressure space
-  pres_regularization.clear();
-//  DoFTools::make_hanging_node_constraints( dof_handler_pressure, pres_regularization ); //??
-
-  // Add the only constraint that we have
-  pres_regularization.add_line(0);
-
-  // close it
-  pres_regularization.close();
-
-  // condense the sparsity pattern
-  pres_regularization.condense( spar_pattern_pressure );
-
-  // Compress the sparsity pattern for the pressure
-  spar_pattern_pressure.compress();
-
-  //// Init the matrices for the pressure
-  pres_Laplace.reinit( spar_pattern_pressure );
-  pres_Mass.reinit( spar_pattern_pressure );
-
-  // Now we assemble the matrices
-  // The Laplace operator is the projection matrix
-  MatrixCreator::create_laplace_matrix( dof_handler_pressure, quadrature_pressure, pres_Laplace );
-  // The pressure mass matrix to do the pressure update
-  MatrixCreator::create_mass_matrix(  dof_handler_pressure, quadrature_pressure, pres_Mass );
-
-  // Finally we condense the Laplace operator
-  pres_regularization.condense( pres_Laplace );
-}
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::init_gradient_operator(){
-  //// Init the sparsity pattern for the gradient operator
-  spar_pattern_pres_vel.reinit( dof_handler_velocity.n_dofs(), dof_handler_pressure.n_dofs(),
-                                dof_handler_velocity.max_couplings_between_dofs() );
-  DoFTools::make_sparsity_pattern( dof_handler_velocity, dof_handler_pressure, spar_pattern_pres_vel );
-  spar_pattern_pres_vel.compress();
-
-  /*
-      To assemble each component of the gradient operator
-      we need to make a loop over all cells and compute the products of
-          velocity * d_#(pressure)
-      where # is the current space component. For this reason we need two cell
-      iterators, one for the velocity and one for the pressure space.
-  */
-  typename DoFHandler<dim>::active_cell_iterator cell_init       = dof_handler_velocity.begin_active(),
-                                                   cell_end        = dof_handler_velocity.end(),
-                                                   cell,
-                                                   bogus_cell_init = dof_handler_pressure.begin_active(),
-                                                   bogus_cell;
-  /*
-    The FEValues extractors.
-    For the velocity we need values,
-    For the pressure we only need gradients.
-  */
-  FEValues<dim> fe_values_velocity( fe_velocity, quadrature_velocity, update_values | update_JxW_values ),
-                fe_values_pressure( fe_pressure, quadrature_velocity, update_gradients );
-
-  // Usual, and useful, abreviations
-  const unsigned int vel_dofs_per_cell  = fe_velocity.dofs_per_cell,
-                     pres_dofs_per_cell = fe_pressure.dofs_per_cell,
-                     n_q_points         = quadrature_velocity.size();
-
-  // The local gradient operator
-  FullMatrix<double> local_grad( vel_dofs_per_cell, pres_dofs_per_cell );
-
-  // Local to global DoF's map
-  std::vector<unsigned int> vel_local_dof_indices( vel_dofs_per_cell ), pres_local_dof_indices( pres_dofs_per_cell );
-
-  for( unsigned int d=0; d<dim; ++d ){
-    // we reinit the given component of the gradient operator
-    pres_Diff[d].reinit( spar_pattern_pres_vel );
-
-    // The usual cell loop
-    for( cell = cell_init, bogus_cell = bogus_cell_init; cell not_eq cell_end; ++cell, ++bogus_cell ){
-      fe_values_velocity.reinit( cell );
-      fe_values_pressure.reinit( bogus_cell );
-
-      cell->get_dof_indices( vel_local_dof_indices );
-      bogus_cell->get_dof_indices( pres_local_dof_indices );
-
-      // local contributions
-      local_grad = 0.;
-      for( unsigned int q=0; q<n_q_points; ++q ){
-        for( unsigned int i=0; i<vel_dofs_per_cell; ++i )
-          for( unsigned int j=0; j<pres_dofs_per_cell; ++j )
-            local_grad( i, j ) += fe_values_velocity.JxW(q)*fe_values_velocity.shape_value( i, q )
-                                  *fe_values_pressure.shape_grad( j, q)[d];
-      }
-
-      // Add it to the global contributions
-      for( unsigned int i=0; i<vel_dofs_per_cell; ++i )
-        for( unsigned int j=0; j<pres_dofs_per_cell; ++j)
-          pres_Diff[d].add( vel_local_dof_indices[i], pres_local_dof_indices[j], local_grad( i, j)  );
-    }
-  }
-}
-
-
-
-// Create Triangulation method
-template<int dim> void Navier_Stokes_Projection<dim>::Create_Triangulation( const unsigned int n_of_refines ){
-  // A disk
-  GridGenerator::hyper_ball( triangulation );
-  static const HyperBallBoundary<dim> boundary;
-  triangulation.set_boundary( 0, boundary );
-
-/*  // Our domain is a unit square.
-  GridGenerator::hyper_cube(triangulation);*/
-  triangulation.refine_global( n_of_refines );
-  std::cout<<" Number of active cells: "<<triangulation.n_active_cells()<<std::endl;
-
-  // Having created the mesh, we create the mesh-dependent data
-
-  //// Distribute the degrees of freedom and renumber them
-  dof_handler_velocity.distribute_dofs( fe_velocity );
-  DoFRenumbering::boost::Cuthill_McKee( dof_handler_velocity );
-  dof_handler_pressure.distribute_dofs( fe_pressure );
-  DoFRenumbering::boost::Cuthill_McKee( dof_handler_pressure );
-
-  // Initialize the matrices for the velocity
-  init_velocity_matrices();
-
-  // Initialize the matrices for the pressure
-  init_pressure_matrices();
-
-  // Initialize the gradient operator
-  init_gradient_operator();
-
-
-  //// Set the correct size for the various vectors that we need
-  pres_n.reinit( dof_handler_pressure.n_dofs() );
-  pres_n_minus_1.reinit( dof_handler_pressure.n_dofs() );
-
-  phi_n.reinit( dof_handler_pressure.n_dofs() );
-  phi_n_minus_1.reinit( dof_handler_pressure.n_dofs() );
-
-  pres_tmp.reinit( dof_handler_pressure.n_dofs() );
-
-  for(unsigned int d=0; d<dim; ++d){
-    u_n[d].reinit( dof_handler_velocity.n_dofs() );
-    u_n_minus_1[d].reinit( dof_handler_velocity.n_dofs() );
-    u_star[d].reinit( dof_handler_velocity.n_dofs() );
-    force[d].reinit( dof_handler_velocity.n_dofs() );
-  }
-  v_tmp.reinit( dof_handler_velocity.n_dofs() );
-
-  std::cout<<" dim( X_h ) = "<<( dof_handler_velocity.n_dofs()*dim )<<std::endl
-           <<" dim( M_h ) = "<<dof_handler_pressure.n_dofs()<<std::endl;
-}
-
-
-
-// Initialization of the time dependent data
-template<int dim> void Navier_Stokes_Projection<dim>::Initialize(){
-
-  // The constant part of the matrix for the velocity is 1.5/dt*Mass + 1/Re*Delta_h
-  vel_Laplace_plus_Mass = 0.;
-  vel_Laplace_plus_Mass.add( 1./Re, vel_Laplace );
-  vel_Laplace_plus_Mass.add( 1.5/dt, vel_Mass );
-
-  // Now we start the vectors and load them with the initial data
-  ////pressure
-  Pressure<dim> pres( t_0 );
-  VectorTools::interpolate( dof_handler_pressure, pres, pres_n_minus_1 );
-  pres.advance_time( dt );
-  VectorTools::interpolate( dof_handler_pressure, pres, pres_n );
-
-  ////phi
-  phi_n = 0.;
-  phi_n_minus_1 = 0.;
-
-  ////velocity
-  Velocity<dim> v;
-  for(unsigned int d=0; d<dim; ++d){
-    v.set_time( t_0 );
-    v.set_component(d);
-
-    VectorTools::interpolate( dof_handler_velocity, v, u_n_minus_1[d] );
-
-    v.advance_time( dt );
-    VectorTools::interpolate( dof_handler_velocity, v, u_n[d] );
-  }
-  plot_solution(1);
-}
-
-
-
-/*
-    A time marching step, the output messages and name of the functions are self
-    explanatory, so no comments here
-*/
-template<int dim> void Navier_Stokes_Projection<dim>::run( const bool verbose, const unsigned int n_of_plots ){
-  // We set the number of steps
-  unsigned int n_steps = ( T - t_0 )/dt;
-  // Set the source term to the correct time
-  rhs.set_time( 2.*dt );
-  vel_exact.set_time( 2.*dt );
-  // do as many steps as necessary
-  timeval init_time, end_time;
-
-  /// debug
-  show = false;
-  Solve_time = 0.;
-  ///-----
-
-  gettimeofday( &init_time, 0 );
-  for( unsigned int n = 2; n<=n_steps; ++n ){
-
-    /// debug
-    if( n == n_steps )
-      show = true;
-    ///----
-
-    if( verbose )
-      std::cout<<" Step = "<<n<<" Time = "<<(n*dt)<<std::endl;
-    if( verbose )
-      std::cout<<"  Interpolating the velocity "<<std::endl;
-    interpolate_velocity();
-    if( verbose )
-      std::cout<<"  Diffusion Step"<<std::endl;
-    if( ( n%vel_update_prec==0 ) and ( verbose ) )
-      std::cout<<"   With reinitialization of the preconditioner"<<std::endl;
-    diffusion_step( (n%vel_update_prec == 0 ) or ( n == 2) );
-    if( verbose )
-      std::cout<<"  Projection Step"<<std::endl;
-    projection_step( ( n == 2 ) );
-    if( verbose )
-      std::cout<<"  Updating the Pressure"<<std::endl;
-    update_pressure( ( n == 2 ) );
-    if( n%n_of_plots ==0 ){
-      if( verbose )
-        std::cout<<" Plotting Solution"<<std::endl;
-      plot_solution(n);
-    }
-    //  advance in time
-    rhs.advance_time(dt);
-    vel_exact.advance_time(dt);
-  }
-  gettimeofday( &end_time, 0 );
-  double run_time = double( end_time.tv_sec - init_time.tv_sec ) + 1e-6*double( end_time.tv_usec - init_time.tv_usec );
-  run_time /= n_steps;
-  run_time /= double( dim*dof_handler_velocity.n_dofs() + dof_handler_pressure.n_dofs() );
-//   if( verbose )
-    std::cout<<" Run time = "<<run_time<<std::endl;
-}
-
-
-
-// Simple procedure that interpolates the velocity
-template<int dim> void Navier_Stokes_Projection<dim>::interpolate_velocity(){
-  for( unsigned int d=0; d<dim; ++d ){
-    u_star[d] = 0.;
-    u_star[d].add( 2., u_n[d], -1., u_n_minus_1[d] );
-  }
-}
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::assemble_advection_term( const unsigned int d ){
-  // Cell iterator
-  typename DoFHandler<dim>::active_cell_iterator cell     = dof_handler_velocity.begin_active(),
-                                                   cell_end = dof_handler_velocity.end();
-
-  /*
-    The FEValues extractors.
-    We need values because of the terms involving u^*
-    also we need the gradient to assemble the advection part
-  */
-  FEValues<dim> fe_values( fe_velocity, quadrature_velocity, update_values | update_JxW_values | update_gradients );
-
-  // Usual, and useful, abreviations
-  const unsigned int dofs_per_cell = fe_velocity.dofs_per_cell,
-                     n_q_points    = quadrature_velocity.size();
-
-  /*
-    We need to compute the values of several finite element functions (and some gradients)
-    at quadrature points. For this we use all these arrays.
-    The names are self-explanatory
-  */
-  std::vector<double> u_old_local( n_q_points );
-  std::vector< Point<dim> > u_star_local( n_q_points );
-  std::vector< Tensor<1,dim> > grad_u_star(n_q_points);
-
-  // Local to global DoF's map
-  std::vector<unsigned int> local_dof_indices( dofs_per_cell );
-
-  // Local contribution of the advection operator
-  Vector<double> local_rhs( dofs_per_cell );
-
-  // loop over cells
-  for( ; cell not_eq cell_end; ++cell ){
-    // reinit the stuff
-    fe_values.reinit( cell );
-    cell->get_dof_indices( local_dof_indices );
-
-    // create the local u^*
-    for( unsigned int s=0; s<dim; ++s ){
-      fe_values.get_function_values( u_star[s], u_old_local );
-      for( unsigned int q=0; q<n_q_points; ++q )
-        u_star_local[q](s) = u_old_local[q];
-    }
-
-    fe_values.get_function_gradients( u_star[d], grad_u_star );
-
-    local_rhs = 0.;
-    for( unsigned int q =0; q<n_q_points; ++q )
-      for( unsigned int i=0; i<dofs_per_cell; ++i )
-        // The local contribution of the advection term
-        local_rhs(i) -= ( ( u_star_local[q]*grad_u_star[q] )*fe_values.shape_value( i, q ) )*fe_values.JxW(q) ;
-
-    // Add local contributions to the global matrix
-    for( unsigned int i=0; i<dofs_per_cell; ++i )
-      force[d]( local_dof_indices[i] ) += local_rhs(i);
-  }
-}
-
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::threaded_solve( const unsigned int d){
-
-  // Create the solver that we are going to use
-  SolverControl solver_control( vel_max_its, vel_eps*force[d].l2_norm() );
-  SolverCG<> cg( solver_control );
-  // Solve
-  cg.solve( vel_it_matrix[d], u_n[d], force[d], prec_velocity[d] );
-}
-
-
-
-// A diffusion step
-template<int dim> void Navier_Stokes_Projection<dim>::diffusion_step( const bool reinit_prec ){
-  // Used to solve in parallel the components of the velocity
-  Threads::ThreadGroup<> threads;
-
-  /// debug
-  static double assembly_time = 0.;
-  timeval init_time_assembly, end_time_assembly, init_time_solve, end_time_solve;
-  gettimeofday( &init_time_assembly, 0 );
-  ///----
-
-  /*
-    Define the pressure interpolant. It is defined by
-        p^* = p^{n} + 4/3 phi^{n} - 1/3 phi^{n-1}
-  */
-  pres_tmp = pres_n;
-  pres_tmp.add(4./3., phi_n, -1./3., phi_n_minus_1);
-  pres_tmp *= -1.;
-
-  for( unsigned int d=0; d<dim; ++d ){
-    // First we need to construct the right hand sides
-    rhs.set_component(d);
-    //// compute the source term
-    VectorTools::create_right_hand_side( dof_handler_velocity, quadrature_velocity, rhs, force[d] );
-
-    // We assemble the advection term
-    threads += Threads::spawn( *this, &Navier_Stokes_Projection<dim>::assemble_advection_term )(d);
-    //assemble_advection_term(d);
-
-    /*
-      Compute the value of the contribution of the old velocities.
-      Recall that BDF2 reads
-      ( 0.5/dt)*( 3u^{n+1} - 4u^{n} + u^{n-1} )
-      the terms involving u^{n} and u^{n-1} are known, so they go into the
-      right hand side
-    */
-    v_tmp = 0.;
-    v_tmp.add( 2./dt,u_n[d],-.5/dt,u_n_minus_1[d] );
-    vel_Mass.vmult_add( force[d], v_tmp );
-
-    // And add the contribution of the pressure part
-    pres_Diff[d].vmult_add( force[d], pres_tmp );
-
-    // Copy the old value
-    u_n_minus_1[d] = u_n[d];
-
-
-    // Create the matrix
-    vel_it_matrix[d].copy_from( vel_Laplace_plus_Mass );
-    //// Create the boundary values
-    vel_exact.set_component(d);
-    VectorTools::interpolate_boundary_values( dof_handler_velocity, 0, vel_exact, boundary_values );
-
-    // apply BCs
-    MatrixTools::apply_boundary_values( boundary_values, vel_it_matrix[d], u_n[d], force[d] );
-  }
-  threads.join_all();
-
-  ///  debug
-  gettimeofday( &end_time_assembly, 0 );
-  assembly_time += double( end_time_assembly.tv_sec - init_time_assembly.tv_sec )
-                            + 1e-6*double( end_time_assembly.tv_usec - init_time_assembly.tv_usec );
-  gettimeofday( &init_time_solve, 0 );
-  ///-----
-
-  for(unsigned int d=0; d<dim; ++d ){
-    // Check if the preconditioner needs reinitialization
-  //  if( reinit_prec )
-      prec_velocity[d].initialize( vel_it_matrix[d]/*, SparseILU<double>::AdditionalData( vel_diag_strength, vel_off_diagonals )*/ );
-
-    //prec_velocity[d].factorize( vel_it_matrix[d] );
-//    threaded_solve(d);
-
-//     And solve
-//    threads += Threads::spawn( *this, &Navier_Stokes_Projection<dim>::threaded_solve )(d);
-
-//     prec_velocity[d].solve( force[d] );
-//     u_n[d] = force[d];
-    prec_velocity[d].vmult( force[d], u_n[d] );
-  }
-//   threads.join_all();
-
-
-  ///debug
-  gettimeofday( &end_time_solve, 0 );
-  Solve_time += double( end_time_solve.tv_sec - init_time_solve.tv_sec )
-                            + 1e-6*double( end_time_solve.tv_usec - init_time_solve.tv_usec );
-  if(show)
-    std::cout<<" Assembly time = "<<assembly_time<<std::endl
-             <<" Solve time = "<<Solve_time<<std::endl;
-  ///----
-}
-
-
-
-// A projection step
-template<int dim> void Navier_Stokes_Projection<dim>::projection_step( const bool reinit_prec ){
-
-  /// debug
-  timeval init_time_assemble, init_time_solve, end_time_assemble, end_time_solve, init_time_factorize, end_time_factorize;
-  static double assemble_time=0., solve_time=0., factorize_time = 0.;
-  ///----
-
-
-  // Check if the preconditioner needs reinitialization
-  if( reinit_prec ){
-
-    /// debug
-    gettimeofday( &init_time_factorize, 0 );
-    ///-----
-// prec_pressure.initialize( pres_Laplace, SparseILU<double>::AdditionalData( proj_diag_strength, proj_off_diagonals ) );
-    prec_pressure.initialize( pres_Laplace );
-
-    /// debug
-    gettimeofday( &end_time_factorize, 0 );
-    factorize_time += double( end_time_factorize.tv_sec - init_time_factorize.tv_sec )
-                  + 1e-6*double( end_time_factorize.tv_usec - init_time_factorize.tv_usec );
-    ///--------
-  }
-
-  /// debug
-  gettimeofday( &init_time_assemble, 0 );
-  ///-----
-
-  // To make things faster we assembled the discrete gradient operator
-  // so we just need to multiply by its transpose
-  pres_tmp = 0.;
-  for( unsigned d=0; d<dim; ++d )
-    pres_Diff[d].Tvmult_add( pres_tmp, u_n[d] );
-
-  /// debug
-  gettimeofday( &end_time_assemble, 0 );
-  assemble_time += double( end_time_assemble.tv_sec - init_time_assemble.tv_sec )
-                 + 1e-6*double( end_time_assemble.tv_usec - init_time_assemble.tv_usec );
-  ///-----
-
-  // Up to this point, pres_tmp has the term < u_n, grad q >
-  // we will assemble this only once, so we
-  // Copy the old value of phi
-  phi_n_minus_1 = phi_n;
-
-  phi_n = pres_tmp;
-  phi_n *= 1.5/dt;
-
-  /// debug
-  gettimeofday( &init_time_solve, 0 );
-  ///-------
-
-  // condense because of the constraint
-  pres_regularization.condense( phi_n );
-
-
-  prec_pressure.solve( phi_n );
-  // The method that we are going to use
-//   SolverControl solver_control( proj_max_its, proj_eps*pres_tmp.l2_norm() );
-//   SolverCG<> cg( solver_control );
-
-
-  // And solve
-//   cg.solve( pres_Laplace, phi_n, pres_tmp, prec_pressure );
-
-  // After solving we need to distribute the constrained stuff
-  pres_regularization.distribute( phi_n );
-
-  /// debug
-  gettimeofday( &end_time_solve, 0 );
-  solve_time += double( end_time_solve.tv_sec - init_time_solve.tv_sec )
-                 + 1e-6*double( end_time_solve.tv_usec - init_time_solve.tv_usec );
-  if(show)
-    std::cout<<" DIV assembly time = "<<assemble_time<<std::endl
-             <<" LAP solve time = "<<solve_time<<std::endl
-             <<" LAP factorize time = "<<factorize_time<<std::endl;
-  ///------
-}
-
-
-
-// Pressure update step
-template<int dim> void Navier_Stokes_Projection<dim>::update_pressure( const bool reinit_prec ){
-
-  /// debug
-  static double solve_time = 0.;
-  timeval init_time_solve, end_time_solve;
-  ///---------
-
-  // Copy the old pressure
-  pres_n_minus_1 = pres_n;
-  switch( type ){
-    case METHOD_STANDARD:
-      /*
-        Simple standard case
-            p^{n+1} = p^n + \phi^{n+1}
-      */
-      pres_n += phi_n;
-      break;
-    case METHOD_ROTATIONAL:
-    {
-      /*
-        Rotational form
-            p^{n+1} = p^n + \phi^{n+1} - 1/Re div u^{n+1}
-      */
-
-      /// debug
-      gettimeofday( &init_time_solve, 0 );
-      ///------
-
-      // The preconditioner for the mass matrix problem
-      if( reinit_prec )
-        prec_mass.initialize( pres_Mass/*, SparseILU<double>::AdditionalData( pres_diag_strength, pres_off_diagonals )*/ );
-
-      // The solver for  the mass matrix
-//       SolverControl solver_control( pres_max_its, pres_eps*pres_tmp.l2_norm() );
-//       SolverCG<> cg( solver_control );
-//
-//       // Solve the mass matrix problem
-//       cg.solve( pres_Mass, pres_n, pres_tmp, prec_mass );
-
-      pres_n = pres_tmp;
-      prec_mass.solve( pres_n );
-
-      pres_n.sadd(1./Re, 1., pres_n_minus_1, 1., phi_n );
-
-      /// debug
-      gettimeofday( &end_time_solve, 0 );
-      solve_time += double( end_time_solve.tv_sec - init_time_solve.tv_sec )
-                 + 1e-6*double( end_time_solve.tv_usec - init_time_solve.tv_usec );
-      ///----------
-
-      break;
-    }
-    default:
-      Assert( false, ExcNotImplemented() );
-  };
-
-  /// debug
-  if(show)
-    std::cout<<" PRES_UPD solve time = "<<solve_time<<std::endl;
-  ///---
-
-}
-
-
-
-// Post Processing
-template<int dim> void Navier_Stokes_Projection<dim>::Post_Process(){
-  double tmp, vel_err_L2=0., vel_err_H1=0., pres_err_L2;
-
-  Vector<double> differences( triangulation.n_active_cells() );
-
-  Velocity<dim> u_exact(T);
-  for( unsigned int d=0; d<dim; ++d ){
-    u_exact.set_component(d);
-
-    // L2 error for this component of the velocity
-    differences = 0.;
-    VectorTools::integrate_difference( dof_handler_velocity, u_n[d], u_exact, differences, quadrature_velocity, VectorTools::L2_norm );
-    tmp = differences.l2_norm();
-    vel_err_L2 += tmp*tmp;
-
-    // H1 error for this component of the velocity
-    differences = 0.;
-    VectorTools::integrate_difference( dof_handler_velocity, u_n[d], u_exact, differences, quadrature_velocity, VectorTools::H1_seminorm );
-    tmp = differences.l2_norm();
-    vel_err_H1 += tmp*tmp;
-  }
-  vel_err_L2 = std::sqrt( vel_err_L2 );
-  vel_err_H1 = std::sqrt( vel_err_H1 );
-
-
-  /*
-    The actual pressure needs to have mean value zero.
-    We compute the mean value and subtract it from the computed pressure
-  */
-  double pres_mean_value = VectorTools::compute_mean_value( dof_handler_pressure, quadrature_pressure, pres_n, 0 );
-  for( unsigned int i=0; i<pres_n.size(); ++i )
-    pres_n(i) -= pres_mean_value;
-
-  Pressure<dim> pres_exact(T);
-
-  // L2 error of the pressure
-  differences = 0.;
-  VectorTools::integrate_difference( dof_handler_pressure, pres_n, pres_exact, differences, quadrature_pressure, VectorTools::L2_norm );
-  pres_err_L2 = differences.l2_norm();
-
-  // Add the obtained results to the convergence table
-  convergence_table.add_value( "dt"     , dt );
-  convergence_table.add_value( "u_L2"   , vel_err_L2 );
-  convergence_table.add_value( "u_H1"   , vel_err_H1 );
-  convergence_table.add_value( "pres_L2", pres_err_L2 );
-
-  convergence_table.set_precision( "dt"     , 5 );
-  convergence_table.set_precision( "u_L2"   , 5 );
-  convergence_table.set_precision( "u_H1"   , 5 );
-  convergence_table.set_precision( "pres_L2", 5 );
-
-  convergence_table.set_scientific( "u_L2"   , true );
-  convergence_table.set_scientific( "u_H1"   , true );
-  convergence_table.set_scientific( "pres_L2", true );
-
-  // And show it
-  convergence_table.write_text(std::cout);
-}
-
-
-
-// Output of the solution
-template<int dim> void Navier_Stokes_Projection<dim>::plot_solution( const unsigned int step ){
-  // This only works in 2d
-  Assert( dim==2, ExcNotImplemented() );
-
-  // We need to assemble the vorticy
-  FEValues<dim> fe_values( fe_velocity, quadrature_velocity, update_values | update_gradients | update_JxW_values );
-  const unsigned int dofs_per_cell = fe_velocity.n_dofs_per_cell(),
-                     n_q_points    = quadrature_velocity.size();
-  std::vector< Tensor<1,dim> > grad_vel_1( n_q_points ), grad_vel_2( n_q_points );
-  std::vector< unsigned int> local_dof_indices( dofs_per_cell );
-  Vector<double> local_rhs( dofs_per_cell );
-  double vorticity;
-
-  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler_velocity.begin_active(),
-                                                   cend = dof_handler_velocity.end();
-
-  force[0] = 0.;
-  //We start the usual loop
-  for( ; cell not_eq cend; ++cell ){
-    // reinit all the needed stuff
-    local_rhs = 0.;
-    fe_values.reinit( cell );
-    cell->get_dof_indices( local_dof_indices );
-
-    // get the gradients of each function
-    fe_values.get_function_gradients( u_n[0], grad_vel_1 );
-    fe_values.get_function_gradients( u_n[1], grad_vel_2 );
-
-    // usual loop over quad points and local dofs
-    for( unsigned int q=0; q<n_q_points; ++q ){
-      vorticity = grad_vel_2[q][1] - grad_vel_1[q][0];
-      for( unsigned int i=0; i<dofs_per_cell; ++i )
-        local_rhs(i) += fe_values.JxW(q)*vorticity*fe_values.shape_value( i, q );
-    }
-
-    // And copy it to the local one
-    for( unsigned int i=0; i<dofs_per_cell; ++i )
-      force[0]( local_dof_indices[i] ) += local_rhs(i);
-  }
-
-  // We solve a mass matrix problem
-  SolverControl solver_control( 500, 1e-6*force[0].l2_norm() );
-  SolverCG<> cg( solver_control );
-  static bool is_prec_initted = false;
-  static SparseILU<double> prec;
-  if( not is_prec_initted ){
-    prec.initialize( vel_Mass, SparseILU<double>::AdditionalData( 1e-5, 70 ) );
-    is_prec_initted = true;
-  }
-  cg.solve( vel_Mass, force[1], force[0], prec );
-
-  // Once we have the vorticity we can output it
-  DataOut<dim> data_out;
-  data_out.attach_dof_handler( dof_handler_velocity );
-  data_out.add_data_vector( force[1], "vorticity" );
-  data_out.build_patches();
-  std::ostringstream filename;
-  filename<<"vorticity"<<step<<".vtk";
-  std::ofstream output( filename.str().c_str() );
-  data_out.write_vtk( output );
-}
-
-
-// explicit template instantiation
-template class Navier_Stokes_Projection<deal_II_dimension>;
diff --git a/deal.II/examples/step-35/source/main.cc b/deal.II/examples/step-35/source/main.cc
deleted file mode 100644 (file)
index 69d9738..0000000
+++ /dev/null
@@ -1,66 +0,0 @@
-#include "../include/NavierStokes.h"
-
-
-#include <iostream>
-#include <sys/time.h>
-
-
-
-int main( int argc, char **argv ){
-  try{
-    Data_Storage data;
-    if( argc<2 ){
-      std::cout<<std::endl<<std::endl
-               <<" Missing parameter file!!"<<std::endl
-               <<" Aborting."<<std::endl<<std::endl;
-      data.print_usage();
-      exit(1);
-    }
-    data.read_data( argv[1] );
-    deallog.depth_console( data.verbose?2:0 );
-    Navier_Stokes_Projection<deal_II_dimension> test( data );
-    test.Create_Triangulation( data.n_of_global_refines );
-    timeval init_time, end_time;
-    gettimeofday( &init_time, 0 );
-    for( double dt = data.initial_dt; dt >= data.final_dt; dt /= data.dt_decrement ){
-      std::cout<<" dt = "<<dt<<std::endl;
-      test.set_dt( dt );
-      test.Initialize();
-      test.run( data.verbose, data.output );
-      test.Post_Process();
-      std::cout<<"====================================="<<std::endl<<std::endl;
-    }
-    gettimeofday( &end_time, 0 );
-    std::cout<<" Global time = "
-             <<( double( end_time.tv_sec - init_time.tv_sec )+1e-6*double( end_time.tv_usec - init_time.tv_usec ) )
-             <<" sec"<<std::endl;
-  }
-  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;
-  }
-  std::cout<<"----------------------------------------------------"
-           <<std::endl
-           <<"Apparently everything went fine!"
-           <<std::endl
-           <<"Don't forget to brush your teeth :-)"
-           <<std::endl<<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.