From: bangerth Date: Mon, 12 May 2008 21:26:26 +0000 (+0000) Subject: Go a bit further. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b3eb671e99b07b89299ceba10340d44cb9cf7f88;p=dealii-svn.git Go a bit further. git-svn-id: https://svn.dealii.org/trunk@16084 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-33/doc/intro.dox b/deal.II/examples/step-33/doc/intro.dox index f0e0cc8812..319b6095f1 100644 --- a/deal.II/examples/step-33/doc/intro.dox +++ b/deal.II/examples/step-33/doc/intro.dox @@ -251,13 +251,22 @@ tends not to slow the simulation to a halt. This, of course, is purely a heuris if the author's advisor heard about it, the author would likely be exiled forever from the finite element error estimation community. -

Input Deck

+

Input deck, initial and boundary conditions

We use an input file deck to drive the simulation. In this way, we can alter the boundary conditions and other important properties of the simulation without having to recompile. For more information on the format, look at the results section, where we describe an example input file in more detail. +In previous example programs, we have usually hard-coded the initial +and boundary conditions. In this program, we instead use the +expression parser class FunctionParser so that we can specify a +generic expression in the input file and have it parsed at run time — +this way, we can change initial conditions without the need to +recompile the program. Consequently, no classes named +InitialConditions or BoundaryConditions will be declared in the +program below. + diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index a49ae86710..02a1d38935 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -234,13 +234,19 @@ struct EulerEquations } - // On the boundaries of the domain and across hanging nodes we use - // a numerical flux function to enforce boundary conditions. This routine - // is the basic Lax-Friedrich's flux with a stabilization parameter - // $\alpha$. + // On the boundaries of the + // domain and across hanging + // nodes we use a numerical flux + // function to enforce boundary + // conditions. This routine is + // the basic Lax-Friedrich's flux + // with a stabilization parameter + // $\alpha$. It's form has also + // been given already in the + // introduction: template static - void numerical_normal_flux(const Point &normal, + void numerical_normal_flux(const Point &normal, const std::vector &Wplus, const std::vector &Wminus, const double alpha, @@ -268,127 +274,6 @@ template const double EulerEquations::gas_gamma = 1.4; - - - - // @sect3{Initial and side condition parsing} - // For the initial condition we use the expression parser function - // object. -template -class InitialCondition : public FunctionParser -{ - public: - InitialCondition (); - - // This function should be called after parsing, but before using - // the object. It formalizes the expressions and initializes the - // function parser with the appropriate expressions. - void Init(); - - // During parsing we call this function as the initial condition - // for one of the $\mathbf{w}$ variables is encountered. - - void set_ic(int _row, std::string &expr) { - expressions[_row] = expr; - } - - virtual void vector_value_list (const std::vector > &points, - std::vector > &value_list) const; - private: - std::vector expressions; -}; - -template -InitialCondition::InitialCondition () : - FunctionParser (EulerEquations::n_components), - expressions(EulerEquations::n_components, "0.0") -{} - - // Here we set up x,y,z as the variables that one should use in the input - // deck to describe their initial condition. -template -void InitialCondition::Init() { - std::map constants; - constants["M_PI"] = M_PI; - std::string variables = (dim == 2 ? "x,y" : "x,y,z"); - - FunctionParser::initialize(variables, expressions, constants); - -} - -template -void InitialCondition::vector_value_list (const std::vector > &points, - std::vector > &value_list) const -{ - const unsigned int n_points = points.size(); - - Assert (value_list.size() == n_points, - ExcDimensionMismatch (value_list.size(), n_points)); - - for (unsigned int p=0; p::vector_value (points[p], - value_list[p]); -} - - // As above, we use the expression function parser for boundary conditions. -template -class SideCondition : public FunctionParser -{ - public: - SideCondition (int ncomp); - ~SideCondition (); - - // As above. - void Init(); - // As above. - void set_coeff_row(int _row_n, std::string &expr); - - virtual void vector_value_list (const std::vector > &points, - std::vector > &value_list) const; - private: - std::vector expressions; -}; - -template -SideCondition::SideCondition (int ncomp) : - FunctionParser (ncomp), - expressions(ncomp, "0.0") -{ -} -template -void SideCondition::set_coeff_row (int _row_n, std::string &expr) -{ - expressions[_row_n] = expr; -} - -template -void SideCondition::Init() { - std::map constants; - constants["M_PI"] = M_PI; - std::string variables = (dim == 2 ? "x,y" : "x,y,z"); - - FunctionParser::initialize(variables, expressions, constants); - -} - -template -SideCondition::~SideCondition () -{ -} - -template -void SideCondition::vector_value_list (const std::vector > &points, - std::vector > &value_list) const -{ - const unsigned int n_points = points.size(); - - Assert (value_list.size() == n_points, - ExcDimensionMismatch (value_list.size(), n_points)); - - for (unsigned int p=0; p::vector_value (points[p], - value_list[p]); -} // @sect3{Conservation Law class} // Here we define a Conservation Law class that helps group // operations and data for our Euler equations into a manageable @@ -470,7 +355,7 @@ class ConsLaw ParameterHandler prm; // Name of the mesh to read in. string mesh; - InitialCondition ic; + FunctionParser initial_conditions; // Enums for the various supported boundary conditions. typedef enum {INFLOW_BC = 1, OUTFLOW_BC=2, NO_PENETRATION_BC=3, PRESSURE_BC=4} bc_type; @@ -578,10 +463,8 @@ void ConsLaw::add_boundary(unsigned int bd, template void ConsLaw::initialize() { VectorTools::interpolate(dof_handler, - ic, solution); - VectorTools::interpolate(dof_handler, - ic, nlsolution); - + initial_conditions, solution); + nlsolution = solution; } // @sect3{Assembly} @@ -1228,6 +1111,7 @@ ConsLaw::ConsLaw () T(0), dT(0.05), TF(10), + initial_conditions (EulerEquations::n_components), is_stationary(false), Map(NULL), Matrix(NULL), @@ -1877,12 +1761,15 @@ void ConsLaw::load_parameters(const char *infile){ // Define a parser for every boundary, though it may be // unused. - SideCondition *sd = new SideCondition(EulerEquations::n_components); + FunctionParser *sd = new FunctionParser(EulerEquations::n_components); + + std::vector expressions(EulerEquations::n_components, "0.0"); + char bd[512]; std::sprintf(bd, "boundary_%d", b); prm.enter_subsection(bd); - const std::string &nopen = prm.get("no penetration"); + const std::string nopen = prm.get("no penetration"); // Determine how each component is handled. for (unsigned int di = 0; di < EulerEquations::n_components; di++) { @@ -1896,29 +1783,35 @@ void ConsLaw::load_parameters(const char *infile){ flags[di] = NO_PENETRATION_BC; } else if (btype == "inflow") { flags[di] = INFLOW_BC; - sd->set_coeff_row(di, var_value); + expressions[di] = var_value; } else if (btype == "pressure") { flags[di] = PRESSURE_BC; - sd->set_coeff_row(di, var_value); + expressions[di] = var_value; } } prm.leave_subsection(); // Add the boundary condition to the law. - sd->Init(); + sd->initialize (FunctionParser::default_variable_names(), + expressions, + std::map()); add_boundary(b, flags, sd); } // Initial conditions. prm.enter_subsection("initial condition"); - for (unsigned int di = 0; di < EulerEquations::n_components; di++) { - char var[512]; - - std::sprintf(var, "w_%d value", di); - std::string var_value = prm.get(var); - ic.set_ic(di, var_value); + { + std::vector expressions; + for (unsigned int di = 0; di < EulerEquations::n_components; di++) + { + char var[512]; + std::sprintf(var, "w_%d value", di); + expressions.push_back (prm.get(var)); + } + initial_conditions.initialize (FunctionParser::default_variable_names(), + expressions, + std::map()); } - ic.Init(); prm.leave_subsection(); // The linear solver. @@ -2138,13 +2031,18 @@ void ConsLaw::run () // The following ``main'' function is // similar to previous examples and - // need not to be commented on. + // need not to be commented on. Note + // that the program aborts if no + // input file name is given on the + // command line. int main (int argc, char *argv[]) { - if (argc != 2) { - std::cout << "Usage:" << argv[0] << " infile" << std::endl; - std::exit(1); - } + if (argc != 2) + { + std::cout << "Usage:" << argv[0] << " infile" << std::endl; + std::exit(1); + } + try { ConsLaw<2> cons;