From: bangerth Date: Fri, 18 Sep 2009 00:55:58 +0000 (+0000) Subject: Reindent a bunch of stuff. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=79ed3eb03b61b57b7904234f28226ac7b99aefd0;p=dealii-svn.git Reindent a bunch of stuff. git-svn-id: https://svn.dealii.org/trunk@19473 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-35/step-35.cc b/deal.II/examples/step-35/step-35.cc index e59fee8310..ec48e880ac 100644 --- a/deal.II/examples/step-35/step-35.cc +++ b/deal.II/examples/step-35/step-35.cc @@ -61,31 +61,43 @@ #include #include -// Finally we import all the deal.II names to the global namespace + // Finally we import all the deal.II + // names to the global namespace using namespace dealii; -// @sect3{Run time parameters} -// -// Since our method has several parameters that can be fine-tuned we put them -// into an external file, so that they can be determined at run-time. -// -// This includes, in particular, the formulation of the equation for the auxiliary variable $\phi$, for which we declare an enum. -// Next, we declare a class that is going to read and store all the parameters that our program -// needs to run. -namespace RunTimeParameters{ - enum Method_Formulation{ + // @sect3{Run time parameters} + // + // Since our method has several + // parameters that can be fine-tuned + // we put them into an external file, + // so that they can be determined at + // run-time. + // + // This includes, in particular, the + // formulation of the equation for + // the auxiliary variable $\phi$, for + // which we declare an + // enum. Next, we + // declare a class that is going to + // read and store all the parameters + // that our program needs to run. +namespace RunTimeParameters +{ + enum MethodFormulation + { METHOD_STANDARD, METHOD_ROTATIONAL }; - class Data_Storage{ + class Data_Storage + { public: Data_Storage(); ~Data_Storage(); - void read_data( const char *filename ); - Method_Formulation form; + void read_data (const char *filename); + MethodFormulation form; double initial_time, final_time, Reynolds; @@ -104,253 +116,306 @@ namespace RunTimeParameters{ ParameterHandler prm; }; -// In the constructor of this class we declare all the parameters. The -// details of how this works have been discussed elsewhere, for example in -// step-19 and step-29. - 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. " ); - prm.enter_subsection( "Physical data" ); + // In the constructor of this class + // we declare all the + // parameters. The details of how + // this works have been discussed + // elsewhere, for example in + // step-19 and step-29. + 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. "); + 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.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(); - prm.enter_subsection( "Time step data" ); + prm.enter_subsection ("Time step data"); { - prm.declare_entry( "dt", "5e-4", Patterns::Double( 0. ), " The time step size. " ); + prm.declare_entry ("dt", "5e-4", + Patterns::Double (0.), + " The time step size. "); } prm.leave_subsection(); - prm.enter_subsection( "Space discretization" ); + prm.enter_subsection ("Space discretization"); { - prm.declare_entry( "n_of_refines", "0", Patterns::Integer( 0, 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.declare_entry ("n_of_refines", "0", + Patterns::Integer (0, 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(); - prm.enter_subsection( "Data solve velocity" ); + 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.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(); - prm.declare_entry( "verbose", "true", Patterns::Bool(), - " This indicates whether the output of the solution process should be verbose. " ); + prm.declare_entry ("verbose", "true", + Patterns::Bool(), + " This indicates whether the output of the solution process should be verbose. "); - prm.declare_entry( "output", "1", Patterns::Integer(1), - " This indicates between how many time steps we print the solution. " ); + prm.declare_entry ("output", "1", + Patterns::Integer(1), + " This indicates between how many time steps we print the solution. "); } - Data_Storage::~Data_Storage(){} + Data_Storage::~Data_Storage() + {} - void Data_Storage::read_data( const char *filename ){ - std::ifstream file( filename ); - AssertThrow (file, ExcFileNotOpen( filename )); + void Data_Storage::read_data (const char *filename) + { + std::ifstream file (filename); + AssertThrow (file, ExcFileNotOpen (filename)); - prm.read_input( file ); + prm.read_input (file); - if( prm.get( "Method_Form" ) == std::string( "rotational" ) ) + if (prm.get ("Method_Form") == std::string ("rotational")) form = METHOD_ROTATIONAL; else form = METHOD_STANDARD; - prm.enter_subsection( "Physical data" ); + prm.enter_subsection ("Physical data"); { - initial_time = prm.get_double( "initial_time" ); - final_time = prm.get_double( "final_time" ); - Reynolds = prm.get_double( "Reynolds" ); + initial_time = prm.get_double ("initial_time"); + final_time = prm.get_double ("final_time"); + Reynolds = prm.get_double ("Reynolds"); } prm.leave_subsection(); - prm.enter_subsection( "Time step data" ); + prm.enter_subsection ("Time step data"); { - dt = prm.get_double( "dt" ); + dt = prm.get_double ("dt"); } prm.leave_subsection(); - prm.enter_subsection( "Space discretization" ); + prm.enter_subsection ("Space discretization"); { - n_of_global_refines = prm.get_integer( "n_of_refines" ); - pressure_degree = prm.get_integer( "pressure_fe_degree" ); + n_of_global_refines = prm.get_integer ("n_of_refines"); + pressure_degree = prm.get_integer ("pressure_fe_degree"); } prm.leave_subsection(); - prm.enter_subsection( "Data solve velocity" ); + 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" ); + 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(); - verbose = prm.get_bool( "verbose" ); + verbose = prm.get_bool ("verbose"); - output = prm.get_integer( "output" ); + output = prm.get_integer ("output"); } } -// @sect3{The Equation Data} -// In the next namespace, we declare the initial and boundary conditions: -namespace EquationData{ - // As we have chosen a completely decoupled formulation, we will not take advantage of deal.II's capabilities - // to handle vector valued problems. We do, however, - // want to use an interface for the equation data that is - // somehow dimension independent. To be able to do that, our functions should be able to know - // on which spatial component we are currently working, and we should be able to have a - // common interface to do that. The following class is an attempt in that direction. - template - class MultiComponentFunction: public Function{ + // @sect3{The Equation Data} + + // In the next namespace, we declare + // the initial and boundary + // conditions: +namespace EquationData +{ + // As we have chosen a completely + // decoupled formulation, we will + // not take advantage of deal.II's + // capabilities to handle vector + // valued problems. We do, however, + // want to use an interface for the + // equation data that is somehow + // dimension independent. To be + // able to do that, our functions + // should be able to know on which + // spatial component we are + // currently working, and we should + // be able to have a common + // interface to do that. The + // following class is an attempt in + // that direction. + template + class MultiComponentFunction: public Function + { public: - MultiComponentFunction( const double initial_time = 0. ); - void set_component( const unsigned int d ); + MultiComponentFunction (const double initial_time = 0.); + void set_component (const unsigned int d); protected: unsigned int comp; }; - template + template MultiComponentFunction:: - MultiComponentFunction( const double initial_time ) + MultiComponentFunction (const double initial_time) : - Function( 1, initial_time ), comp(0) + Function (1, initial_time), comp(0) {} - template - void MultiComponentFunction::set_component(const unsigned int d ) + template + void MultiComponentFunction::set_component(const unsigned int d) { - Assert( d + // With this class defined, we + // declare classes that describe + // the boundary conditions for + // velocity and pressure: + template class Velocity : public MultiComponentFunction { public: - Velocity( const double initial_time = 0.0 ); + Velocity (const double initial_time = 0.0); - virtual double value( const Point &p, - const unsigned int component = 0 ) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - virtual void value_list( const std::vector< Point > &points, + virtual void value_list (const std::vector< Point > &points, std::vector &values, - const unsigned int component = 0 ) const; + const unsigned int component = 0) const; }; - template - Velocity::Velocity( const double initial_time ) + template + Velocity::Velocity (const double initial_time) : - MultiComponentFunction( initial_time ) + MultiComponentFunction (initial_time) {} - template - void Velocity::value_list( const std::vector > &points, - std::vector &values, - const unsigned int ) const + template + void Velocity::value_list (const std::vector > &points, + std::vector &values, + const unsigned int) const { const unsigned int n_points = points.size(); - Assert( values.size() == n_points, - ExcDimensionMismatch( values.size(), n_points ) ); + Assert (values.size() == n_points, + ExcDimensionMismatch (values.size(), n_points)); for (unsigned int i=0; i::value( points[i] ); + values[i] = Velocity::value (points[i]); } - template - double Velocity::value( const Point &p, - const unsigned int ) const + template + double Velocity::value (const Point &p, + const unsigned int) const { - if( this->comp == 0 ) + if (this->comp == 0) { const double Um = 1.5; const double H = 4.1; - return 4.*Um*p(1)*( H - p(1) )/(H*H); + return 4.*Um*p(1)* (H - p(1))/(H*H); } else return 0; } - template + + + template class Pressure: public Function { public: - Pressure( const double initial_time = 0.0 ); + Pressure (const double initial_time = 0.0); - virtual double value( const Point &p, - const unsigned int component = 0 ) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - virtual void value_list( const std::vector< Point > &points, + virtual void value_list (const std::vector< Point > &points, std::vector &values, - const unsigned int component = 0 ) const; + const unsigned int component = 0) const; }; - template - Pressure::Pressure( const double initial_time ) + template + Pressure::Pressure (const double initial_time) : - Function( 1, initial_time ) + Function (1, initial_time) {} - template - double Pressure::value( const Point &p, - const unsigned int ) const + template + double Pressure::value (const Point &p, + const unsigned int) const { return 0.; } - template - void Pressure::value_list( const std::vector > &points, + template + void Pressure::value_list (const std::vector > &points, std::vector &values, - const unsigned int ) const + const unsigned int) const { const unsigned int n_points = points.size(); - Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) ); + Assert (values.size() == n_points, ExcDimensionMismatch (values.size(), n_points)); for (unsigned int i=0; i::value( points[i] ); + values[i] = Pressure::value (points[i]); } } -// @sect3{The Navier_Stokes_Projection class} -// Now for the main class of the program. It implements the various versions of the projection -// method for Navier-Stokes equations. -// The names for all the methods and attributes should be self-explanatory. -template class Navier_Stokes_Projection{ + // @sect3{The NavierStokesProjection class} + + // Now for the main class of the + // program. It implements the various + // versions of the projection method + // for Navier-Stokes equations. The + // names for all the methods and + // attributes should be + // self-explanatory. +template +class NavierStokesProjection +{ public: - Navier_Stokes_Projection( const RunTimeParameters::Data_Storage &data ); - ~Navier_Stokes_Projection(); - void run( const bool verbose = false, const unsigned int n_of_plots = 10 ); + NavierStokesProjection (const RunTimeParameters::Data_Storage &data); + ~NavierStokesProjection(); + void run (const bool verbose = false, + const unsigned int n_of_plots = 10); protected: - RunTimeParameters::Method_Formulation type; + RunTimeParameters::MethodFormulation type; unsigned int deg; double dt; @@ -366,7 +431,7 @@ template class Navier_Stokes_Projection{ FE_Q fe_velocity, fe_pressure; QGauss quadrature_pressure, quadrature_velocity; - SparsityPattern spar_pattern_velocity, spar_pattern_pressure, spar_pattern_pres_vel; + SparsityPattern sparsity_pattern_velocity, sparsity_pattern_pressure, sparsity_pattern_pres_vel; SparseMatrix vel_Laplace_plus_Mass, vel_it_matrix[dim], vel_Mass, vel_Laplace, vel_Advection, pres_Laplace, pres_Mass, pres_Diff[dim]; @@ -380,16 +445,16 @@ template class Navier_Stokes_Projection{ SparseILU prec_velocity[dim]; SparseDirectUMFPACK prec_mass, prec_pressure, prec_vel_mass; - DeclException2( ExcInvalidTimeStep, double, double, <<" The time step "< class Navier_Stokes_Projection{ typedef std_cxx1x::tuple< typename DoFHandler::active_cell_iterator, typename DoFHandler::active_cell_iterator > IteratorTuple; - typedef parallel::internal::SynchronousIterators SIterators; - struct InitGradPerTaskData{ - unsigned int d, vel_dpc, pres_dpc; - FullMatrix local_grad; - std::vector vel_local_dof_indices, pres_local_dof_indices; - InitGradPerTaskData( const unsigned int dd, const unsigned int vdpc, const unsigned int pdpc ): - d(dd), vel_dpc( vdpc ), pres_dpc( pdpc ), - local_grad( vdpc, pdpc ), vel_local_dof_indices( vdpc ), - pres_local_dof_indices( pdpc ){} - }; - struct InitGradScratchData{ - unsigned int nqp; - FEValues fe_val_vel, fe_val_pres; - InitGradScratchData( const FE_Q &fe_v, const FE_Q &fe_p, const QGauss &quad, - const UpdateFlags flags_v, const UpdateFlags flags_p ) : - nqp( quad.size() ), fe_val_vel( fe_v, quad, flags_v ), - fe_val_pres( fe_p, quad, flags_p ){} - InitGradScratchData( const InitGradScratchData &data ): nqp( data.nqp ), - fe_val_vel( data.fe_val_vel.get_fe(), data.fe_val_vel.get_quadrature(), - data.fe_val_vel.get_update_flags() ), - fe_val_pres( data.fe_val_pres.get_fe(), data.fe_val_pres.get_quadrature(), - data.fe_val_pres.get_update_flags() ) {} - }; - void assemble_one_cell_of_gradient( const SIterators &SI, InitGradScratchData &scratch, - InitGradPerTaskData &data ); - void copy_gradient_local_to_global( const InitGradPerTaskData &data ); - - void assemble_advection_term(); - struct AdvectionPerTaskData{ - FullMatrix local_advection; - std::vector local_dof_indices; - AdvectionPerTaskData( const unsigned int dpc ): local_advection( dpc, dpc ), local_dof_indices( dpc ) {} - }; - struct AdvectionScratchData{ - unsigned int nqp, dpc; - std::vector< Point > u_star_local; - std::vector< Tensor<1,dim> > grad_u_star; - std::vector u_star_tmp; - FEValues fe_val; - AdvectionScratchData( const FE_Q &fe, const QGauss &quad, const UpdateFlags flags ): - nqp( quad.size() ), dpc( fe.dofs_per_cell ), - u_star_local( nqp ), grad_u_star( nqp ), u_star_tmp( nqp ), - fe_val( fe, quad, flags ){} - AdvectionScratchData( const AdvectionScratchData &data ): nqp( data.nqp ), dpc( data.dpc ), - u_star_local( nqp ), grad_u_star( nqp ), u_star_tmp( nqp ), - fe_val( data.fe_val.get_fe(), data.fe_val.get_quadrature(), - data.fe_val.get_update_flags() ) {} - }; +typedef parallel::internal::SynchronousIterators SIterators; +struct InitGradPerTaskData +{ + unsigned int d, vel_dpc, pres_dpc; + FullMatrix local_grad; + std::vector vel_local_dof_indices, pres_local_dof_indices; + InitGradPerTaskData (const unsigned int dd, const unsigned int vdpc, const unsigned int pdpc): + d(dd), vel_dpc (vdpc), pres_dpc (pdpc), + local_grad (vdpc, pdpc), vel_local_dof_indices (vdpc), + pres_local_dof_indices (pdpc) + {} +}; +struct InitGradScratchData +{ + unsigned int nqp; + FEValues fe_val_vel, fe_val_pres; + InitGradScratchData (const FE_Q &fe_v, const FE_Q &fe_p, const QGauss &quad, + const UpdateFlags flags_v, const UpdateFlags flags_p) : + nqp (quad.size()), fe_val_vel (fe_v, quad, flags_v), + fe_val_pres (fe_p, quad, flags_p) + {} + InitGradScratchData (const InitGradScratchData &data): nqp (data.nqp), + fe_val_vel (data.fe_val_vel.get_fe(), data.fe_val_vel.get_quadrature(), + data.fe_val_vel.get_update_flags()), + fe_val_pres (data.fe_val_pres.get_fe(), data.fe_val_pres.get_quadrature(), + data.fe_val_pres.get_update_flags()) + {} +}; +void assemble_one_cell_of_gradient (const SIterators &SI, InitGradScratchData &scratch, + InitGradPerTaskData &data); +void copy_gradient_local_to_global (const InitGradPerTaskData &data); + +void assemble_advection_term(); +struct AdvectionPerTaskData +{ + FullMatrix local_advection; + std::vector local_dof_indices; + AdvectionPerTaskData (const unsigned int dpc): local_advection (dpc, dpc), local_dof_indices (dpc) + {} +}; +struct AdvectionScratchData +{ + unsigned int nqp, dpc; + std::vector< Point > u_star_local; + std::vector< Tensor<1,dim> > grad_u_star; + std::vector u_star_tmp; + FEValues fe_val; + AdvectionScratchData (const FE_Q &fe, const QGauss &quad, const UpdateFlags flags): + nqp (quad.size()), dpc (fe.dofs_per_cell), + u_star_local (nqp), grad_u_star (nqp), u_star_tmp (nqp), + fe_val (fe, quad, flags) + {} + AdvectionScratchData (const AdvectionScratchData &data) + : + nqp (data.nqp), dpc (data.dpc), + u_star_local (nqp), grad_u_star (nqp), u_star_tmp (nqp), + fe_val (data.fe_val.get_fe(), data.fe_val.get_quadrature(), + data.fe_val.get_update_flags()) + {} +}; - void assemble_one_cell_of_advection( const typename DoFHandler::active_cell_iterator &cell, - AdvectionScratchData &scratch, AdvectionPerTaskData &data ); - void copy_advection_local_to_global( const AdvectionPerTaskData &data ); - inline void diffusion_component_solve( const unsigned int d ); +void assemble_one_cell_of_advection (const typename DoFHandler::active_cell_iterator &cell, + AdvectionScratchData &scratch, AdvectionPerTaskData &data); +void copy_advection_local_to_global (const AdvectionPerTaskData &data); +inline void diffusion_component_solve (const unsigned int d); - inline void plot_solution( const unsigned int step ); - inline void assemble_vorticity( const bool reinit_prec ); +inline void plot_solution (const unsigned int step); +inline void assemble_vorticity (const bool reinit_prec); }; -template Navier_Stokes_Projection::~Navier_Stokes_Projection(){ +template +NavierStokesProjection::~NavierStokesProjection() +{ dof_handler_velocity.clear(); dof_handler_pressure.clear(); } -// @sect4{ Navier_Stokes_Projection::Navier_Stokes_Projection } -// In the constructor, we just read all the data from the Data_Storage -// object that is passed as an argument, verify that the data we read is reasonable -// and, finally, create the triangulation and load the initial data. -template Navier_Stokes_Projection::Navier_Stokes_Projection( - const RunTimeParameters::Data_Storage &data ): - type( data.form ), deg( data.pressure_degree ), dt( data.dt ), t_0( data.initial_time ), - T( data.final_time ), Re( data.Reynolds ), 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) + // @sect4{ NavierStokesProjection::NavierStokesProjection } + + // In the constructor, we just read + // all the data from the + // Data_Storage object + // that is passed as an argument, + // verify that the data we read is + // reasonable and, finally, create + // the triangulation and load the + // initial data. +template +NavierStokesProjection::NavierStokesProjection(const RunTimeParameters::Data_Storage &data) + : + type (data.form), deg (data.pressure_degree), dt (data.dt), t_0 (data.initial_time), + T (data.final_time), Re (data.Reynolds), 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) { if(deg < 1) - std::cout<<" WARNING: The chosen pair of finite element spaces is not stable."< .5*T ) ), ExcInvalidTimeStep( dt, .5*T ) ); + AssertThrow (not ( (dt <= 0.) or (dt > .5*T)), ExcInvalidTimeStep (dt, .5*T)); - Create_Triangulation( data.n_of_global_refines ); + Create_Triangulation (data.n_of_global_refines); Initialize(); } -// @sect4{ Navier_Stokes_Projection::Create_Triangulation } -// The method that creates the triangulation and refines it the needed number of times. -// After creating the triangulation, it creates the mesh dependent data, i.e. it distributes -// degrees of freedom and renumbers them, and initializes the matrices and vectors -// that we will use. -template void Navier_Stokes_Projection::Create_Triangulation( const unsigned int n_of_refines ){ + // @sect4{ NavierStokesProjection::Create_Triangulation } + + // The method that creates the + // triangulation and refines it the + // needed number of times. After + // creating the triangulation, it + // creates the mesh dependent data, + // i.e. it distributes degrees of + // freedom and renumbers them, and + // initializes the matrices and + // vectors that we will use. +template +void NavierStokesProjection::Create_Triangulation (const unsigned int n_of_refines) +{ GridIn grid_in; - grid_in.attach_triangulation( triangulation ); + grid_in.attach_triangulation (triangulation); std::string filename = "nsbench2.inp"; - std::ifstream file( filename.c_str() ); - Assert( file, ExcFileNotOpen( filename.c_str() ) ); - grid_in.read_ucd( file ); + std::ifstream file (filename.c_str()); + Assert (file, ExcFileNotOpen (filename.c_str())); + grid_in.read_ucd (file); file.close(); - std::cout<<" Number of refines = "< void Navier_Stokes_Projection::interpolate_velocity(){ - for( unsigned int d=0; d +void NavierStokesProjection::interpolate_velocity() +{ + for (unsigned int d=0; dNavier_Stokes_Projection::diffusion_step} -// The implementation of a diffusion step. -template void Navier_Stokes_Projection::diffusion_step( const bool reinit_prec ){ - pres_tmp.equ( -1., pres_n, -4./3., phi_n, 1./3., phi_n_minus_1 ); + // @sect4{NavierStokesProjection::diffusion_step} + + // The implementation of a diffusion + // step. +template +void NavierStokesProjection::diffusion_step (const bool reinit_prec) +{ + pres_tmp.equ (-1., pres_n, -4./3., phi_n, 1./3., phi_n_minus_1); assemble_advection_term(); - for( unsigned int d=0; d::const_iterator boundaries = boundary_indicators.begin(), - b_end = boundary_indicators.end(); - boundary_values.clear(); - for( ; boundaries not_eq b_end; ++boundaries ){ - switch( *boundaries ){ - case 1: - VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, - ZeroFunction(), boundary_values ); - break; - case 2: - VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, - vel_exact, boundary_values ); - break; - case 3: - if( d not_eq 0 ) - VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, - ZeroFunction(), boundary_values ); - break; - case 4: - VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, - ZeroFunction(), boundary_values ); - break; - default: - Assert( false, ExcNotImplemented() ); - } + for (unsigned int d=0; d::const_iterator boundaries = boundary_indicators.begin(), + b_end = boundary_indicators.end(); + boundary_values.clear(); + for (; boundaries not_eq b_end; ++boundaries) + { + switch (*boundaries) + { + case 1: + VectorTools::interpolate_boundary_values (dof_handler_velocity, *boundaries, + ZeroFunction(), boundary_values); + break; + case 2: + VectorTools::interpolate_boundary_values (dof_handler_velocity, *boundaries, + vel_exact, boundary_values); + break; + case 3: + if (d not_eq 0) + VectorTools::interpolate_boundary_values (dof_handler_velocity, *boundaries, + ZeroFunction(), boundary_values); + break; + case 4: + VectorTools::interpolate_boundary_values (dof_handler_velocity, *boundaries, + ZeroFunction(), boundary_values); + break; + default: + Assert (false, ExcNotImplemented()); + } + } + MatrixTools::apply_boundary_values (boundary_values, vel_it_matrix[d], u_n[d], force[d]); } - MatrixTools::apply_boundary_values( boundary_values, vel_it_matrix[d], u_n[d], force[d] ); - } Threads::TaskGroup tasks; - for(unsigned int d=0; d::AdditionalData( vel_diag_strength, vel_off_diagonals ) ); - tasks += Threads::new_task( &Navier_Stokes_Projection::diffusion_component_solve, *this, d ); - } + for(unsigned int d=0; d::AdditionalData (vel_diag_strength, vel_off_diagonals)); + tasks += Threads::new_task (&NavierStokesProjection::diffusion_component_solve, *this, d); + } tasks.join_all(); } -template void Navier_Stokes_Projection::diffusion_component_solve( const unsigned int d ){ - SolverControl solver_control( vel_max_its, vel_eps*force[d].l2_norm() ); - SolverGMRES<> gmres( solver_control, SolverGMRES<>::AdditionalData( vel_Krylov_size ) ); - gmres.solve( vel_it_matrix[d], u_n[d], force[d], prec_velocity[d] ); +template +void NavierStokesProjection::diffusion_component_solve (const unsigned int d) +{ + SolverControl solver_control (vel_max_its, vel_eps*force[d].l2_norm()); + SolverGMRES<> gmres (solver_control, SolverGMRES<>::AdditionalData (vel_Krylov_size)); + gmres.solve (vel_it_matrix[d], u_n[d], force[d], prec_velocity[d]); } -// @sect4{ The Navier_Stokes_Projection::assemble_advection_term method and related} -template void Navier_Stokes_Projection::assemble_advection_term(){ + // @sect4{ The NavierStokesProjection::assemble_advection_term method and related} +template +void NavierStokesProjection::assemble_advection_term() +{ vel_Advection = 0.; - AdvectionPerTaskData data( fe_velocity.dofs_per_cell ); - AdvectionScratchData scratch( fe_velocity, quadrature_velocity, - update_values | update_JxW_values | update_gradients ); - WorkStream::run( dof_handler_velocity.begin_active(), dof_handler_velocity.end(), *this, - &Navier_Stokes_Projection::assemble_one_cell_of_advection, - &Navier_Stokes_Projection::copy_advection_local_to_global, scratch, data); + AdvectionPerTaskData data (fe_velocity.dofs_per_cell); + AdvectionScratchData scratch (fe_velocity, quadrature_velocity, + update_values | update_JxW_values | update_gradients); + WorkStream::run (dof_handler_velocity.begin_active(), dof_handler_velocity.end(), *this, + &NavierStokesProjection::assemble_one_cell_of_advection, + &NavierStokesProjection::copy_advection_local_to_global, scratch, data); } -template void Navier_Stokes_Projection::assemble_one_cell_of_advection( +template +void NavierStokesProjection::assemble_one_cell_of_advection( const typename DoFHandler::active_cell_iterator &cell, - AdvectionScratchData &scratch, AdvectionPerTaskData &data ){ + AdvectionScratchData &scratch, AdvectionPerTaskData &data) +{ scratch.fe_val.reinit(cell); - cell->get_dof_indices( data.local_dof_indices ); - for( unsigned int d=0; dget_dof_indices (data.local_dof_indices); + for (unsigned int d=0; d void Navier_Stokes_Projection::copy_advection_local_to_global( - const AdvectionPerTaskData &data ){ - for( unsigned int i=0; i +void NavierStokesProjection::copy_advection_local_to_global( + const AdvectionPerTaskData &data) +{ + for (unsigned int i=0; iNavier_Stokes_Projection::projection_step} -// This implements the projection step. -template void Navier_Stokes_Projection::projection_step( const bool reinit_prec ){ - if( reinit_prec ) - prec_pressure.initialize( pres_Laplace ); + // @sect4{NavierStokesProjection::projection_step} + + // This implements the projection step: +template +void NavierStokesProjection::projection_step (const bool reinit_prec) +{ + if (reinit_prec) + prec_pressure.initialize (pres_Laplace); pres_tmp = 0.; - for( unsigned d=0; dNavier_Stokes_Projection::update_pressure } -// This is the pressure update step of the projection method. It implements the -// standard formulation of the method, that is -// @f[ -// p^{n+1} = p^n + \phi^{n+1}, -// @f] -// or the rotational form, which is -// @f[ -// p^{n+1} = p^n + \phi^{n+1} - \frac{1}{Re} \nabla\cdot u^{n+1}. -// @f] -template void Navier_Stokes_Projection::update_pressure( const bool reinit_prec ){ + // @sect4{ NavierStokesProjection::update_pressure } + + // This is the pressure update step + // of the projection method. It + // implements the standard + // formulation of the method, that is + // @f[ + // p^{n+1} = p^n + \phi^{n+1}, + // @f] + // or the rotational form, which is + // @f[ + // p^{n+1} = p^n + \phi^{n+1} - \fraC{1}{Re} \nabla\cdot u^{n+1}. + // @f] +template +void NavierStokesProjection::update_pressure (const bool reinit_prec) +{ pres_n_minus_1 = pres_n; - switch( type ){ - case RunTimeParameters::METHOD_STANDARD: - pres_n += phi_n; - break; - case RunTimeParameters::METHOD_ROTATIONAL: - if( reinit_prec ) - prec_mass.initialize( pres_Mass ); - pres_n = pres_tmp; - prec_mass.solve( pres_n ); - pres_n.sadd(1./Re, 1., pres_n_minus_1, 1., phi_n ); - break; - default: - Assert( false, ExcNotImplemented() ); - }; + switch (type) + { + case RunTimeParameters::METHOD_STANDARD: + pres_n += phi_n; + break; + case RunTimeParameters::METHOD_ROTATIONAL: + if (reinit_prec) + prec_mass.initialize (pres_Mass); + pres_n = pres_tmp; + prec_mass.solve (pres_n); + pres_n.sadd(1./Re, 1., pres_n_minus_1, 1., phi_n); + break; + default: + Assert (false, ExcNotImplemented()); + }; } -// @sect4{ Navier_Stokes_Projection::plot_solution } -// This method plots the current solution. It is an adaptation of -// step-31 and so I will not elaborate on it. -// There is one small detail here. It is often interested to see -// the vorticity of the flow. But, since we are using it here only for plotting -// purposes, we are not going to compute it at every time step, but only when -// we are going to plot it. -template void Navier_Stokes_Projection::plot_solution( const unsigned int step ){ - assemble_vorticity( ( step == 1 ) ); - const FESystem joint_fe( fe_velocity, dim, fe_pressure, 1, fe_velocity, 1 ); - DoFHandler joint_dof_handler( triangulation ); - joint_dof_handler.distribute_dofs( joint_fe ); - Assert( joint_dof_handler.n_dofs() == (dim + 1)*dof_handler_velocity.n_dofs() + dof_handler_pressure.n_dofs(), - ExcInternalError() ); - static Vector joint_solution( joint_dof_handler.n_dofs() ); - std::vector loc_joint_dof_indices( joint_fe.dofs_per_cell ), - loc_vel_dof_indices( fe_velocity.dofs_per_cell ), - loc_pres_dof_indices( fe_pressure.dofs_per_cell ); + // @sect4{ NavierStokesProjection::plot_solution } + + // This method plots the current + // solution. It is an adaptation of + // step-31 and so I will not + // elaborate on it. There is one + // small detail here. It is often + // interested to see the vorticity of + // the flow. But, since we are using + // it here only for plotting + // purposes, we are not going to + // compute it at every time step, but + // only when we are going to plot it. +template +void NavierStokesProjection::plot_solution (const unsigned int step) +{ + assemble_vorticity ( (step == 1)); + const FESystem joint_fe (fe_velocity, dim, fe_pressure, 1, fe_velocity, 1); + DoFHandler joint_dof_handler (triangulation); + joint_dof_handler.distribute_dofs (joint_fe); + Assert (joint_dof_handler.n_dofs() == (dim + 1)*dof_handler_velocity.n_dofs() + dof_handler_pressure.n_dofs(), + ExcInternalError()); + static Vector joint_solution (joint_dof_handler.n_dofs()); + std::vector loc_joint_dof_indices (joint_fe.dofs_per_cell), + loc_vel_dof_indices (fe_velocity.dofs_per_cell), + loc_pres_dof_indices (fe_pressure.dofs_per_cell); typename DoFHandler::active_cell_iterator joint_cell = joint_dof_handler.begin_active(), joint_endc = joint_dof_handler.end(), vel_cell = dof_handler_velocity.begin_active(), pres_cell = dof_handler_pressure.begin_active(); - for( ; joint_cell not_eq joint_endc; ++joint_cell, ++vel_cell, ++pres_cell ){ - joint_cell->get_dof_indices( loc_joint_dof_indices ); - vel_cell->get_dof_indices( loc_vel_dof_indices ), - pres_cell->get_dof_indices( loc_pres_dof_indices ); - for( unsigned int i=0; i joint_solution_names( dim, "v" ); - joint_solution_names.push_back( "p" ); - joint_solution_names.push_back( "rot_u" ); + for (; joint_cell not_eq joint_endc; ++joint_cell, ++vel_cell, ++pres_cell) + { + joint_cell->get_dof_indices (loc_joint_dof_indices); + vel_cell->get_dof_indices (loc_vel_dof_indices), + pres_cell->get_dof_indices (loc_pres_dof_indices); + for (unsigned int i=0; i joint_solution_names (dim, "v"); + joint_solution_names.push_back ("p"); + joint_solution_names.push_back ("rot_u"); DataOut data_out; data_out.attach_dof_handler (joint_dof_handler); std::vector< DataComponentInterpretation::DataComponentInterpretation > - component_interpretation( dim+2, DataComponentInterpretation::component_is_part_of_vector ); + component_interpretation (dim+2, DataComponentInterpretation::component_is_part_of_vector); component_interpretation[dim] = DataComponentInterpretation::component_is_scalar; component_interpretation[dim+1] = DataComponentInterpretation::component_is_scalar; - data_out.add_data_vector( joint_solution, joint_solution_names, DataOut::type_dof_data, - component_interpretation ); - data_out.build_patches( deg + 1 ); + data_out.add_data_vector (joint_solution, joint_solution_names, DataOut::type_dof_data, + component_interpretation); + data_out.build_patches (deg + 1); std::ostringstream filename; - filename<<"solution-"< void Navier_Stokes_Projection::assemble_vorticity( const bool reinit_prec ){ - Assert( dim == 2, ExcNotImplemented() ); - if( reinit_prec ) - prec_vel_mass.initialize( vel_Mass ); + // Since this function is supposed to + // be called only when the plot is + // going to be made, which should not + // be every time step, we do not + // parallelize it. Of course, if + // needed, this can be done as in the + // other cases. Moreover, the + // implementation that we have here + // only works for 2d, so we bail if + // that is not the case. +template +void NavierStokesProjection::assemble_vorticity (const bool reinit_prec) +{ + Assert (dim == 2, ExcNotImplemented()); + if (reinit_prec) + prec_vel_mass.initialize (vel_Mass); typename DoFHandler::active_cell_iterator cell = dof_handler_velocity.begin_active(), end = dof_handler_velocity.end(); - FEValues fe_val_vel( fe_velocity, quadrature_velocity, - update_gradients | update_JxW_values | update_values ); + FEValues fe_val_vel (fe_velocity, quadrature_velocity, + update_gradients | update_JxW_values | update_values); const unsigned int dpc = fe_velocity.dofs_per_cell, nqp = quadrature_velocity.size(); - std::vector ldi( dpc ); - Vector loc_rot( dpc ); + std::vector ldi (dpc); + Vector loc_rot (dpc); - std::vector< Tensor<1,dim> > grad_u1( nqp ), grad_u2( nqp ); + std::vector< Tensor<1,dim> > grad_u1 (nqp), grad_u2 (nqp); rot_u = 0.; - for( ; cell not_eq end; ++cell ){ - fe_val_vel.reinit( cell ); - cell->get_dof_indices( ldi ); - fe_val_vel.get_function_gradients( u_n[0], grad_u1 ); - fe_val_vel.get_function_gradients( u_n[1], grad_u2 ); - loc_rot = 0.; - for( unsigned int q=0; qget_dof_indices (ldi); + fe_val_vel.get_function_gradients (u_n[0], grad_u1); + fe_val_vel.get_function_gradients (u_n[1], grad_u2); + loc_rot = 0.; + for (unsigned int q=0; q test( data ); - test.run( data.verbose, data.output ); - } - 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<<"----------------------------------------------------" - < test (data); + test.run (data.verbose, data.output); + } + 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; }