From: bangerth Date: Fri, 11 Sep 2009 03:19:57 +0000 (+0000) Subject: Minor edits up to the beginning of the main class. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=366282855f497917ba8bbaa92e4096f3aeb00937;p=dealii-svn.git Minor edits up to the beginning of the main class. git-svn-id: https://svn.dealii.org/trunk@19438 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 09237354ce..ce79271619 100644 --- a/deal.II/examples/step-35/step-35.cc +++ b/deal.II/examples/step-35/step-35.cc @@ -13,7 +13,7 @@ // @sect3{Include Files} // We start by including all the necessary - // dealii header files and some C++ related + // deal.II header files and some C++ related // ones. Each one of them has been discussed // in previous tutorial programs, so we will // not get into details here. @@ -26,6 +26,14 @@ #include #include +#include +#include +#include +#include +#include +#include +#include + #include #include #include @@ -45,32 +53,25 @@ #include #include -#include -#include -#include -#include -#include -#include -#include - #include #include #include #include -#include #include #include -// Finally we import all the dealii 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 options that can be fine-tuned we decided to group all these into -// an external file, so that these can be determined at run-time.
-// First, the formulation of the method, which we set as a enum. +// @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{ @@ -103,30 +104,37 @@ namespace RunTimeParameters{ ParameterHandler prm; }; -// In the constructor of this class we declare all the parameters. -// The details of how this works have been discussed somewhere else *** -// so let's not elaborate on that +// 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 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.leave_subsection(); prm.enter_subsection( "Time step data" ); + { prm.declare_entry( "dt", "5e-4", Patterns::Double( 0. ), " The time step size. " ); + } prm.leave_subsection(); 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.leave_subsection(); 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. " ); @@ -137,6 +145,7 @@ namespace RunTimeParameters{ " 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(), @@ -146,65 +155,75 @@ namespace RunTimeParameters{ " This indicates between how many time steps we print the solution. " ); } + + Data_Storage::~Data_Storage(){} + + void Data_Storage::read_data( const char *filename ){ std::ifstream file( filename ); - if( not file ) - throw ExcFileNotOpen( filename ); + AssertThrow (file, ExcFileNotOpen( filename )); + prm.read_input( file ); - std::string token = prm.get( "Method_Form" ); - if( token == std::string( "rotational" ) ) + if( prm.get( "Method_Form" ) == std::string( "rotational" ) ) form = METHOD_ROTATIONAL; else form = METHOD_STANDARD; prm.enter_subsection( "Physical data" ); + { initial_time = prm.get_double( "initial_time" ); final_time = prm.get_double( "final_time" ); Reynolds = prm.get_double( "Reynolds" ); token = prm.get( "filename" ); + } prm.leave_subsection(); prm.enter_subsection( "Time step data" ); + { dt = prm.get_double( "dt" ); + } prm.leave_subsection(); 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(); 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(); verbose = prm.get_bool( "verbose" ); output = prm.get_integer( "output" ); - - file.close(); } } // @sect3{The Equation Data} -// Here we declare the initial and boundary conditions. +// In the next namespace, we declare the initial and boundary conditions: namespace EquationData{ - // Because of our implementation, we do not take advantage of the capabilities of the library - // to handle vector valued problems. Whether this is a good or bad idea is another issue. The - // point here is that we want to be able to write an interface for the equation data that is + // 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 indenpendent. To be able to do that, our functions should be able to know // on which space 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{ + template + class MultiComponentFunction: public Function{ public: MultiComponentFunction( const double initial_time = 0. ); void set_component( const unsigned int d ); @@ -212,58 +231,104 @@ namespace EquationData{ unsigned int comp; }; - template MultiComponentFunction::MultiComponentFunction( const double initial_time ): - Function( 1, initial_time ), comp(0){} - - template void MultiComponentFunction::set_component(const unsigned int d ){ + template + MultiComponentFunction:: + MultiComponentFunction( const double initial_time ) + : + Function( 1, initial_time ), comp(0) + {} + + + template + void MultiComponentFunction::set_component(const unsigned int d ) + { Assert( d class Velocity: public MultiComponentFunction{ + + // 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 ); - virtual double value( const Point &p, const unsigned int component = 0 ) const; - virtual void value_list( const std::vector< Point > &points, std::vector &values, - const unsigned int component = 0 ) const; - }; - template Velocity::Velocity( const double initial_time ): - MultiComponentFunction( initial_time ){} + virtual double value( const Point &p, + const unsigned int component = 0 ) const; - template void Velocity::value_list( const std::vector > &points, - std::vector &values, const unsigned int ) const{ + virtual void value_list( const std::vector< Point > &points, + std::vector &values, + const unsigned int component = 0 ) const; + }; + + + template + Velocity::Velocity( const double initial_time ) + : + MultiComponentFunction( initial_time ) + {} + + + 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] ); } - template inline double Velocity::value( const Point &p, const unsigned int ) const{ + + template + inline double Velocity::value( const Point &p, + const unsigned int ) const + { double return_value = 0.; static const double Um = 1.5, H = 4.1; - if( MultiComponentFunction::comp == 0 ) + if( this->comp == 0 ) return_value = 4.*Um*p(1)*( H - p(1) )/(H*H); return return_value; } - template class Pressure: public Function{ + template + class Pressure: public Function + { public: Pressure( const double initial_time = 0.0 ); - virtual double value( const Point &p, const unsigned int component = 0 ) const; - virtual void value_list( const std::vector< Point > &points, std::vector &values, - 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, + std::vector &values, + const unsigned int component = 0 ) const; }; - template Pressure::Pressure( const double initial_time ): Function( 1, initial_time ){} - - template inline double Pressure::value( const Point &p, const unsigned int ) const{ + template + Pressure::Pressure( const double initial_time ) + : + Function( 1, initial_time ) + {} + + + template + inline double Pressure::value( const Point &p, + const unsigned int ) const + { return 0.; } - template void Pressure::value_list( const std::vector > &points, - std::vector &values, - const unsigned int ) const{ + template + void Pressure::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 ) ); for (unsigned int i=0; iNavier_Stokes_Projection class} -// This is the main class of the program. It implements the various avatars of the projection -// methods for Navier-Stokes equations. -// The names for all the methods and attributes are self-explanatory. +// 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{ public: Navier_Stokes_Projection( const RunTimeParameters::Data_Storage &data );