From: bangerth Date: Tue, 23 Mar 2010 07:06:33 +0000 (+0000) Subject: Delete directories that have been supplanted by the single file step-35.cc. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=70b1cfc1cd78d016a46969becbe9fa21e98a2321;p=dealii-svn.git Delete directories that have been supplanted by the single file step-35.cc. git-svn-id: https://svn.dealii.org/trunk@20885 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-35/include/EqData.h b/deal.II/examples/step-35/include/EqData.h deleted file mode 100644 index d528e5d16c..0000000000 --- a/deal.II/examples/step-35/include/EqData.h +++ /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 - - -// We need to define sines and cosines -#include - - - -// This is ugly, we have to remove it -const double PI = std::acos( -1. ); - - - -// The Velocity function -template class Velocity: public Multi_Component_Function { - public: - Velocity(const double initial_time =0.0); - virtual double value(const Point &p, const unsigned int component = 0) const; - virtual Tensor<1,dim> gradient(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 void gradient_list( const std::vector< Point > &points, std::vector< Tensor<1,dim> > &gradients, - const unsigned int component = 0 ) const; -}; - - - -// The Pressure 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 Tensor<1,dim> gradient(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 void gradient_list( const std::vector< Point > &points, std::vector< Tensor<1,dim> > &gradients, - const unsigned int component = 0 ) const; -}; - - - -// The Force function -template class Force: public Multi_Component_Function{ - public: - Force( 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; -}; - -#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 index 4047aab37e..0000000000 --- a/deal.II/examples/step-35/include/FileReader.h +++ /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 - - - -#include -#include - - - -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 index 22ba2939f4..0000000000 --- a/deal.II/examples/step-35/include/MultiComponentFunction.h +++ /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 - - -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 class Multi_Component_Function: public Function { - 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 index 13c8bfc7f3..0000000000 --- a/deal.II/examples/step-35/include/NavierStokes.h +++ /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 -#include - - - -/* - The includes needed to manage multithreading -*/ -#include -#include - - - -/* - Grid management includes. - Triangulation handler - Grid generator & refinement, accessor and iterators -*/ -#include -#include -#include -#include -#include -#include -#include - - - -/* - Degrees of freedom management includes - Handling, accessing, computing useful stuff from and for them and renumbering. -*/ -#include -#include -#include -#include -#include - - - -/* - Finite element management includes - Definition of finite elements and extraction of info from them - and from finite element functions. -*/ -#include -#include -#include - - - -/* - Linear algebra includes - Sparse matrices, vectors, preconditioners and solvers. -*/ -#include -#include -#include -#include -#include -#include -#include - - - -/* - Numerical algorithms includes - Assembly of right hand sides, computation of errors, output of the solution, etc -*/ -#include -#include -#include - - -/* - This class is the core of all the program - It is the implementation of the projection method for the Navier-Stokes equations. -*/ -template 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 rhs; - //// The exact velocity (to apply boundary values) - Velocity vel_exact; - //// The boundary conditions - std::map boundary_values; - - // Finite Element Spaces data - //// The mesh - Triangulation triangulation; - //// The DoF handlers - DoFHandler dof_handler_velocity, dof_handler_pressure; - //// The polynomial spaces - FE_Q fe_velocity, fe_pressure; - //// The quadrature formulae - QGauss 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 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 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 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 "< Force::Force( const double initial_time ): Multi_Component_Function( initial_time ){ -} - - - -template void Force::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; i::value( points[i] ); -} - - - -template inline double Force::value(const Point &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::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 Velocity::Velocity(const double initial_time): Multi_Component_Function( 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 ) ); - for (unsigned int i=0; i::value( points[i] ); -} - - - -template inline double Velocity::value(const Point &p, const unsigned int) const{ - double return_value = std::cos( Function::get_time() ); - switch( Multi_Component_Function::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 inline Tensor<1,dim> Velocity::gradient(const Point &p, const unsigned int) const{ - Tensor<1,dim> return_value; - switch( Multi_Component_Function::component ){ - // [0, -cos(t)] - case 0: - return_value[0] = 0.; - return_value[1] = -std::cos( Function::get_time() ); - break; - case 1: - // [cos(t), 0] - return_value[0] = std::cos( Function::get_time() ); - return_value[1] = 0.; - break; - default: - Assert( false, ExcNotImplemented() ); - }; - return return_value; -} - - - -template void Velocity::gradient_list( const std::vector > &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::gradient( points[i] ); -} - - - -// Pressure function methods -template Pressure::Pressure(const double initial_time): Function(1,initial_time){} - - - -template inline double Pressure::value(const Point &p, const unsigned int) const{ - // sin(x-y+t) - return std::sin( p(0) )*std::sin( p(1) )*std::sin( Function::get_time() ); -} - - - -template inline Tensor<1,dim> Pressure::gradient(const Point &p, const unsigned int) const{ - // [cos(x)*sin(y)*sin(t), sin(x)*cos(y)*sin(t)] - return Point( std::cos( p(0) )*std::sin( p(1) )*std::sin( Function::get_time() ), - std::sin( p(0) )*std::cos( p(1) )*std::sin( Function::get_time() ) ); -} - - - -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; i::value( points[i] ); -} - - - -template inline void Pressure::gradient_list( const std::vector > &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::gradient( points[i] ); -} - - - -// explicit template instantiation -template class Force; -template class Velocity; -template class Pressure; 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 index 2f52b5acf9..0000000000 --- a/deal.II/examples/step-35/source/EqData_square.ccc +++ /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 Force::Force( const double initial_time ): Multi_Component_Function( initial_time ){ -} - - - -template void Force::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; i::value( points[i] ); -} - - - -template inline double Force::value(const Point &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::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 Velocity::Velocity(const double initial_time): Multi_Component_Function( 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 ) ); - for (unsigned int i=0; i::value( points[i] ); -} - - - -template inline double Velocity::value(const Point &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::get_time() ), - return_value = 0.; - switch( Multi_Component_Function::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 inline Tensor<1,dim> Velocity::gradient(const Point &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::get_time() ); - switch( Multi_Component_Function::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 void Velocity::gradient_list( const std::vector > &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::gradient( points[i] ); -} - - - -// Pressure function methods -template Pressure::Pressure(const double initial_time): Function(1,initial_time){} - - - -template inline double Pressure::value(const Point &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::get_time() ); -} - - - -template inline Tensor<1,dim> Pressure::gradient(const Point &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::get_time() ) )* Point( -std::sin( PI*p(0) )*std::sin( PI*p(1) ), - std::cos( PI*p(0) )*std::cos( PI*p(1) ) ) ; -} - - - -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; i::value( points[i] ); -} - - - -template inline void Pressure::gradient_list( const std::vector > &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::gradient( points[i] ); -} - - - -// explicit template instantiation -template class Force; -template class Velocity; -template class Pressure; diff --git a/deal.II/examples/step-35/source/FileReader.cc b/deal.II/examples/step-35/source/FileReader.cc deleted file mode 100644 index ce58a8c36d..0000000000 --- a/deal.II/examples/step-35/source/FileReader.cc +++ /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" )< Multi_Component_Function::Multi_Component_Function( const double initial_time ): - Function( 1, initial_time ), component(0) { -} - - - -// Set Component Function: Check that it is in range and then set it -template void Multi_Component_Function::set_component(const unsigned int d ){ - // Check if the requested component is correct - Assert( d= 2, ExcNotImplemented() ); - component = d; -} - - - -// explicit template instantiation -template class Multi_Component_Function; diff --git a/deal.II/examples/step-35/source/NavierStokes.cc b/deal.II/examples/step-35/source/NavierStokes.cc deleted file mode 100644 index 88961bd709..0000000000 --- a/deal.II/examples/step-35/source/NavierStokes.cc +++ /dev/null @@ -1,792 +0,0 @@ -/* - Implementation of the Navier_Stokes_projection class - by Abner Salgado. -*/ -#include "../include/NavierStokes.h" - - - -#include - -/// debug -bool show; -double Solve_time; -///---- - - -// Constructor -template Navier_Stokes_Projection::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."< .5*T ) ), ExcInvalidTimeStep( dt, .5*T ) ); -} - - - -// Destructor -template Navier_Stokes_Projection::~Navier_Stokes_Projection(){ - dof_handler_velocity.clear(); - dof_handler_pressure.clear(); -} - - - -// Set time step -template void Navier_Stokes_Projection::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 void Navier_Stokes_Projection::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 void Navier_Stokes_Projection::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 void Navier_Stokes_Projection::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::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 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 local_grad( vel_dofs_per_cell, pres_dofs_per_cell ); - - // Local to global DoF's map - std::vector vel_local_dof_indices( vel_dofs_per_cell ), pres_local_dof_indices( pres_dofs_per_cell ); - - for( unsigned int d=0; dget_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 void Navier_Stokes_Projection::Create_Triangulation( const unsigned int n_of_refines ){ - // A disk - GridGenerator::hyper_ball( triangulation ); - static const HyperBallBoundary 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: "< void Navier_Stokes_Projection::plot_solution( const unsigned int step ){ - // This only works in 2d - Assert( dim==2, ExcNotImplemented() ); - - // We need to assemble the vorticy - FEValues 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 local_rhs( dofs_per_cell ); - double vorticity; - - typename DoFHandler::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 cg( solver_control ); - static bool is_prec_initted = false; - static SparseILU prec; - if( not is_prec_initted ){ - prec.initialize( vel_Mass, SparseILU::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 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"<; diff --git a/deal.II/examples/step-35/source/main.cc b/deal.II/examples/step-35/source/main.cc deleted file mode 100644 index 69d97384a5..0000000000 --- a/deal.II/examples/step-35/source/main.cc +++ /dev/null @@ -1,66 +0,0 @@ -#include "../include/NavierStokes.h" - - -#include -#include - - - -int main( int argc, char **argv ){ - try{ - Data_Storage data; - if( argc<2 ){ - std::cout< 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 = "<