From a328f5908fb673d3551162ada8b1870a07f3e5e9 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 20 Aug 2008 12:06:54 +0000 Subject: [PATCH] Sequester everything that has to do with equation data into its own namespace. git-svn-id: https://svn.dealii.org/trunk@16604 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 394 +++++++++++++--------------- 1 file changed, 184 insertions(+), 210 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index df5e563308..abd377d268 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -294,6 +294,162 @@ void PreconditionerTrilinosAmg::vmult (Vector &dst, } + // @sect3{Equation data} + + // Again, the next stage in the program + // is the definition of the equation + // data, that is, the various + // boundary conditions, the right hand + // side and the initial condition (remember + // that we're about to solve a time- + // dependent system). The basic strategy + // for this definition is the same as in + // step-22. Regarding the details, though, + // there are some differences. + + // The first + // thing is that we don't set any boundary + // conditions on the velocity, as is + // explained in the introduction. So + // what is left are two conditions for + // pressure p and temperature + // T. + + // Secondly, we set an initial + // condition for all problem variables, + // i.e., for u, p and T, + // so the function has dim+2 + // components. + // In this case, we choose a very simple + // test case, where everything is zero. + + // @sect4{Boundary values} +namespace EquationData +{ + // define viscosity + const double eta = 1; + const double kappa = 1e-6; + + template + class PressureBoundaryValues : public Function + { + public: + PressureBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + + template + double + PressureBoundaryValues::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } + + + + + + // @sect4{Initial values} + template + class TemperatureInitialValues : public Function + { + public: + TemperatureInitialValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + }; + + + template + double + TemperatureInitialValues::value (const Point &, + const unsigned int) const + { + return 0; + } + + + template + void + TemperatureInitialValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = TemperatureInitialValues::value (p, c); + } + + + + // @sect4{Right hand side} + // + // The last definition of this kind + // is the one for the right hand + // side function. Again, the content + // of the function is very + // basic and zero in most of the + // components, except for a source + // of temperature in some isolated + // regions near the bottom of the + // computational domain, as is explained + // in the problem description in the + // introduction. + template + class TemperatureRightHandSide : public Function + { + public: + TemperatureRightHandSide () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + }; + + + template + double + TemperatureRightHandSide::value (const Point &p, + const unsigned int /*component*/) const + { + static const Point source_centers[3] + = { (dim == 2 ? Point(.3,.1) : Point(.3,.5,.1)), + (dim == 2 ? Point(.45,.1) : Point(.45,.5,.1)), + (dim == 2 ? Point(.75,.1) : Point(.75,.5,.1)) }; + static const double source_radius + = (dim == 2 ? 1./32 : 1./8); + + return ((source_centers[0].distance (p) < source_radius) + || + (source_centers[1].distance (p) < source_radius) + || + (source_centers[2].distance (p) < source_radius) + ? + 1 + : + 0); + } + + + template + void + TemperatureRightHandSide::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = TemperatureRightHandSide::value (p, c); + } +} + + // @sect3{The BoussinesqFlowProblem class template} @@ -371,188 +527,13 @@ class BoussinesqFlowProblem boost::shared_ptr > Mp_preconditioner; bool rebuild_stokes_matrix; - bool rebuild_temperature_matrix; + bool rebuild_temperature_matrices; bool rebuild_stokes_preconditioner; }; - // @sect3{Equation data} - - // Again, the next stage in the program - // is the definition of the equation - // data, that is, the various - // boundary conditions, the right hand - // side and the initial condition (remember - // that we're about to solve a time- - // dependent system). The basic strategy - // for this definition is the same as in - // step-22. Regarding the details, though, - // there are some differences. - - // The first - // thing is that we don't set any boundary - // conditions on the velocity, as is - // explained in the introduction. So - // what is left are two conditions for - // pressure p and temperature - // T. - - // Secondly, we set an initial - // condition for all problem variables, - // i.e., for u, p and T, - // so the function has dim+2 - // components. - // In this case, we choose a very simple - // test case, where everything is zero. - - // @sect4{Boundary values} -template -class PressureBoundaryValues : public Function -{ - public: - PressureBoundaryValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - -template -double -PressureBoundaryValues::value (const Point &/*p*/, - const unsigned int /*component*/) const -{ - return 0; -} - - - -template -class TemperatureBoundaryValues : public Function -{ - public: - TemperatureBoundaryValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - - -template -double -TemperatureBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const -{ -//TODO: leftover from olden times. replace by something sensible once we have -//diffusion in the temperature field - if (p[0] == 0) - return 1; - else - return 0; -} - - - - // @sect4{Initial values} -template -class TemperatureInitialValues : public Function -{ - public: - TemperatureInitialValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; -}; - - -template -double -TemperatureInitialValues::value (const Point &, - const unsigned int) const -{ - return 0; -} - - -template -void -TemperatureInitialValues::vector_value (const Point &p, - Vector &values) const -{ - for (unsigned int c=0; cn_components; ++c) - values(c) = TemperatureInitialValues::value (p, c); -} - - - - // @sect4{Right hand side} - // - // The last definition of this kind - // is the one for the right hand - // side function. Again, the content - // of the function is very - // basic and zero in most of the - // components, except for a source - // of temperature in some isolated - // regions near the bottom of the - // computational domain, as is explained - // in the problem description in the - // introduction. -template -class TemperatureRightHandSide : public Function -{ - public: - TemperatureRightHandSide () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; -}; - - -template -double -TemperatureRightHandSide::value (const Point &p, - const unsigned int /*component*/) const -{ - static const Point source_centers[3] - = { (dim == 2 ? Point(.3,.1) : Point(.3,.5,.1)), - (dim == 2 ? Point(.45,.1) : Point(.45,.5,.1)), - (dim == 2 ? Point(.75,.1) : Point(.75,.5,.1)) }; - static const double source_radius - = (dim == 2 ? 1./32 : 1./8); - - return ((source_centers[0].distance (p) < source_radius) - || - (source_centers[1].distance (p) < source_radius) - || - (source_centers[2].distance (p) < source_radius) - ? - 1 - : - 0); -} - - -template -void -TemperatureRightHandSide::vector_value (const Point &p, - Vector &values) const -{ - for (unsigned int c=0; cn_components; ++c) - values(c) = TemperatureRightHandSide::value (p, c); -} - - - // @sect3{Linear solvers and preconditioners} @@ -823,7 +804,7 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) old_time_step (0), timestep_number (0), rebuild_stokes_matrix (true), - rebuild_temperature_matrix (true), + rebuild_temperature_matrices (true), rebuild_stokes_preconditioner (true) {} @@ -1471,7 +1452,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () // individual blocks (velocity, // pressure, temperature) from // the total FE system. - const PressureBoundaryValues pressure_boundary_values; + const EquationData::PressureBoundaryValues pressure_boundary_values; std::vector boundary_values (n_face_q_points); std::vector old_temperature_values(n_q_points); @@ -1555,12 +1536,11 @@ void BoussinesqFlowProblem::assemble_stokes_system () } } - // define viscosity - const double eta = 1; if (rebuild_stokes_matrix) for (unsigned int i=0; i > &old_old_temperature_hessians, const std::vector > &present_stokes_values, const std::vector &gamma_values, - const double kappa, const double global_u_infty, const double global_T_variation, const double global_Omega_diameter, @@ -1728,7 +1707,7 @@ double compute_viscosity( const double u_grad_T = u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2; - const double kappa_Delta_T = kappa + const double kappa_Delta_T = EquationData::kappa * (trace(old_temperature_hessians[q]) + trace(old_old_temperature_hessians[q])) / 2; @@ -1778,7 +1757,7 @@ double compute_viscosity( template void BoussinesqFlowProblem::assemble_temperature_matrix () { - if (rebuild_temperature_matrix == false) + if (rebuild_temperature_matrices == false) return; temperature_mass_matrix = 0; @@ -1821,8 +1800,6 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () temperature_fe_values.reinit (cell); - const double kappa = 1e-6; - for (unsigned int q=0; q::assemble_temperature_matrix () for (unsigned int i=0; i::assemble_temperature_matrix () temperature_constraints.condense (temperature_mass_matrix); temperature_constraints.condense (temperature_stiffness_matrix); - rebuild_temperature_matrix = false; + rebuild_temperature_matrices = false; } @@ -1924,8 +1905,7 @@ void BoussinesqFlowProblem::assemble_temperature_system () std::vector > old_old_temperature_hessians(n_q_points); - TemperatureBoundaryValues temperature_boundary_values; - TemperatureRightHandSide temperature_right_hand_side; + EquationData::TemperatureRightHandSide temperature_right_hand_side; std::vector gamma_values (n_q_points); std::vector phi_T (dofs_per_cell); @@ -1980,7 +1960,6 @@ void BoussinesqFlowProblem::assemble_temperature_system () // of diffusion (determined // impirically) to keep the // scheme stable - const double kappa = 1e-6; const double nu = compute_viscosity (old_temperature_values, old_old_temperature_values, @@ -1990,7 +1969,7 @@ void BoussinesqFlowProblem::assemble_temperature_system () old_old_temperature_hessians, present_stokes_values, gamma_values, - kappa, global_u_infty, + global_u_infty, global_T_range.second - global_T_range.first, global_Omega_diameter, cell->diameter(), old_time_step); @@ -2269,19 +2248,17 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(max_grid_level); cell != triangulation.end(); ++cell) - if (cell->has_children() == false) - cell->clear_refine_flag (); + cell->clear_refine_flag (); - SolutionTransfer soltrans(temperature_dof_handler); - - triangulation.prepare_coarsening_and_refinement(); - std::vector > x_solution (2); x_solution[0].reinit (temperature_dof_handler.n_dofs()); x_solution[0] = temperature_solution; x_solution[1].reinit (temperature_dof_handler.n_dofs()); x_solution[1] = old_temperature_solution; + SolutionTransfer soltrans(temperature_dof_handler); + + triangulation.prepare_coarsening_and_refinement(); soltrans.prepare_for_coarsening_and_refinement(x_solution); triangulation.execute_coarsening_and_refinement (); @@ -2296,7 +2273,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) old_temperature_solution = tmp[1]; rebuild_stokes_matrix = true; - rebuild_temperature_matrix = true; + rebuild_temperature_matrices = true; rebuild_stokes_preconditioner = true; } @@ -2306,12 +2283,10 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) template double BoussinesqFlowProblem::get_maximal_velocity () const { - QGauss quadrature_formula(stokes_degree+2); - const unsigned int n_q_points - = quadrature_formula.size(); + const QGauss quadrature_formula(stokes_degree+2); + const unsigned int n_q_points = quadrature_formula.size(); - FEValues fe_values (stokes_fe, quadrature_formula, - update_values); + FEValues fe_values (stokes_fe, quadrature_formula, update_values); std::vector > stokes_values(n_q_points, Vector(dim+1)); double max_velocity = 0; @@ -2330,8 +2305,7 @@ double BoussinesqFlowProblem::get_maximal_velocity () const for (unsigned int i=0; i::run () VectorTools::project (temperature_dof_handler, temperature_constraints, QGauss(temperature_degree+2), - TemperatureInitialValues(), + EquationData::TemperatureInitialValues(), old_temperature_solution); timestep_number = 0; -- 2.39.5