From 7967a969822c0f786c7a7b186e5542fdfc2cdcce Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 20 Mar 2013 12:25:21 +0000 Subject: [PATCH] More cleanups. git-svn-id: https://svn.dealii.org/trunk@28953 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-26/step-26.cc | 143 +++++++++++++++++----------- 1 file changed, 87 insertions(+), 56 deletions(-) diff --git a/deal.II/examples/step-26/step-26.cc b/deal.II/examples/step-26/step-26.cc index 368567c199..eb145c4a55 100644 --- a/deal.II/examples/step-26/step-26.cc +++ b/deal.II/examples/step-26/step-26.cc @@ -41,6 +41,8 @@ namespace Step26 { using namespace dealii; + + template class HeatEquation { @@ -50,38 +52,41 @@ namespace Step26 private: void setup_system(); - void solve_u(); + void solve_time_step(); void output_results() const; - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; - ConstraintMatrix constraints; + ConstraintMatrix constraints; - SparsityPattern sparsity_pattern; + SparsityPattern sparsity_pattern; SparseMatrix mass_matrix; SparseMatrix laplace_matrix; SparseMatrix system_matrix; - Vector solution; - Vector old_solution; - Vector system_rhs; + Vector solution; + Vector old_solution; + Vector system_rhs; - double time, time_step; - unsigned int timestep_number; - const double theta; + double time; + double time_step; + unsigned int timestep_number; + + const double theta; }; - //------------------------------------- + template class RightHandSide: public Function { public: - RightHandSide() : - Function(), - period (0.2) + RightHandSide() + : + Function(), + period (0.2) {} virtual double value(const Point &p, @@ -91,6 +96,8 @@ namespace Step26 const double period; }; + + template double RightHandSide::value(const Point &p, const unsigned int component) const @@ -103,51 +110,60 @@ namespace Step26 if ((point_within_period >= 0.0) && (point_within_period <= 0.2)) { - if ((p[0] > 0.5) && (p[1] > -0.5)) - return 1; - else - return 0; + if ((p[0] > 0.5) && (p[1] > -0.5)) + return 1; + else + return 0; } else if ((point_within_period >= 0.5) && (point_within_period <= 0.7)) { - if ((p[0] > -0.5) && (p[1] > 0.5)) - return 1; - else - return 0; + if ((p[0] > -0.5) && (p[1] > 0.5)) + return 1; + else + return 0; } else return 0; } + + template class BoundaryValues: public Function { public: - BoundaryValues() : + BoundaryValues() + : Function() { } - virtual ~BoundaryValues() - { - } + virtual double value(const Point &p, const unsigned int component = 0) const; }; template double BoundaryValues::value(const Point &/*p*/, - const unsigned int component) const + const unsigned int component) const { Assert(component == 0, ExcInternalError()); - return 0; // Zero-Dirichlet Boundary + return 0; } + + template - HeatEquation::HeatEquation() : - fe(1), dof_handler(triangulation), time_step(1. / 500), theta(0.5) + HeatEquation::HeatEquation() + : + fe(1), + dof_handler(triangulation), + time_step(1. / 500), + theta(0.5) { } + + template void HeatEquation::setup_system() { @@ -162,17 +178,22 @@ namespace Step26 std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl << std::endl; - sparsity_pattern.reinit(dof_handler.n_dofs(), dof_handler.n_dofs(), + sparsity_pattern.reinit(dof_handler.n_dofs(), + dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern); + DoFTools::make_sparsity_pattern(dof_handler, + sparsity_pattern); sparsity_pattern.compress(); mass_matrix.reinit(sparsity_pattern); laplace_matrix.reinit(sparsity_pattern); system_matrix.reinit(sparsity_pattern); - MatrixCreator::create_mass_matrix(dof_handler, QGauss(3), mass_matrix); - MatrixCreator::create_laplace_matrix(dof_handler, QGauss(3), + MatrixCreator::create_mass_matrix(dof_handler, + QGauss(fe.degree+1), + mass_matrix); + MatrixCreator::create_laplace_matrix(dof_handler, + QGauss(fe.degree+1), laplace_matrix); solution.reinit(dof_handler.n_dofs()); @@ -182,8 +203,10 @@ namespace Step26 constraints.close(); } + + template - void HeatEquation::solve_u() + void HeatEquation::solve_time_step() { SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm()); SolverCG<> cg(solver_control); @@ -193,10 +216,12 @@ namespace Step26 cg.solve(system_matrix, solution, system_rhs, preconditioner); - std::cout << " u-equation: " << solver_control.last_step() + std::cout << " " << solver_control.last_step() << " CG iterations." << std::endl; } + + template void HeatEquation::output_results() const { @@ -208,33 +233,34 @@ namespace Step26 data_out.build_patches(); const std::string filename = "solution-" - + Utilities::int_to_string(timestep_number, 3) + ".vtk"; + + Utilities::int_to_string(timestep_number, 3) + + ".vtk"; std::ofstream output(filename.c_str()); data_out.write_vtk(output); - - std::cout << " max= " << time << ' ' << solution.linfty_norm() << std::endl; } + + template void HeatEquation::run() { setup_system(); - VectorTools::interpolate(dof_handler, ZeroFunction(), solution); + VectorTools::interpolate(dof_handler, + ZeroFunction(), + old_solution); + solution = old_solution; timestep_number = 0; output_results(); - VectorTools::interpolate(dof_handler, ZeroFunction(), - old_solution); - Vector tmp(solution.size()); Vector forcing_terms(solution.size()); while (time <= 0.5) { - time += time_step; - ++timestep_number; + time += time_step; + ++timestep_number; std::cout << "Time step " << timestep_number << " at t=" << time << std::endl; @@ -246,14 +272,18 @@ namespace Step26 RightHandSide rhs_function; rhs_function.set_time(time); - VectorTools::create_right_hand_side(dof_handler, QGauss(2), - rhs_function, tmp); + VectorTools::create_right_hand_side(dof_handler, + QGauss(fe.degree+1), + rhs_function, + tmp); forcing_terms = tmp; forcing_terms *= time_step * theta; rhs_function.set_time(time - time_step); - VectorTools::create_right_hand_side(dof_handler, QGauss(2), - rhs_function, tmp); + VectorTools::create_right_hand_side(dof_handler, + QGauss(fe.degree+1), + rhs_function, + tmp); forcing_terms.add(time_step * (1 - theta), tmp); @@ -265,18 +295,19 @@ namespace Step26 std::map boundary_values; VectorTools::interpolate_boundary_values(dof_handler, - 0, + 0, boundary_values_function, - boundary_values); + boundary_values); system_matrix.copy_from(mass_matrix); system_matrix.add(theta * time_step, laplace_matrix); MatrixTools::apply_boundary_values(boundary_values, - system_matrix, + system_matrix, solution, - system_rhs); + system_rhs); } - solve_u(); + + solve_time_step(); output_results(); -- 2.39.5