From 8f79a9e0e0aaa47bf382a2f53337cad5078cb357 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 16 Jul 2010 18:19:21 +0000 Subject: [PATCH] Remove an unused variable. git-svn-id: https://svn.dealii.org/trunk@21510 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-21/step-21.cc | 253 ++++++++++++++-------------- 1 file changed, 126 insertions(+), 127 deletions(-) diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc index 75bf4b7af0..0f4c805e38 100644 --- a/deal.II/examples/step-21/step-21.cc +++ b/deal.II/examples/step-21/step-21.cc @@ -3,7 +3,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2006, 2007, 2008, 2009 by the deal.II authors */ +/* Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -70,7 +70,7 @@ using namespace dealii; // @sect3{The TwoPhaseFlowProblem class} - + // This is the main class of the program. It // is close to the one of step-20, but with a // few additional functions: @@ -101,12 +101,12 @@ using namespace dealii; // several of the formulas in the nonlinear // equations. template -class TwoPhaseFlowProblem +class TwoPhaseFlowProblem { public: TwoPhaseFlowProblem (const unsigned int degree); void run (); - + private: void make_grid_and_dofs (); void assemble_system (); @@ -115,9 +115,9 @@ class TwoPhaseFlowProblem void solve (); void project_back_saturation (); void output_results () const; - + const unsigned int degree; - + Triangulation triangulation; FESystem fe; DoFHandler dof_handler; @@ -126,11 +126,11 @@ class TwoPhaseFlowProblem BlockSparseMatrix system_matrix; const unsigned int n_refinement_steps; - + double time_step; unsigned int timestep_number; - double viscosity; - + double viscosity; + BlockVector solution; BlockVector old_solution; BlockVector system_rhs; @@ -146,11 +146,11 @@ class TwoPhaseFlowProblem // is fully equipped to deal with anything // else, if this is desired: template -class PressureRightHandSide : public Function +class PressureRightHandSide : public Function { public: PressureRightHandSide () : Function(1) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; }; @@ -160,7 +160,7 @@ class PressureRightHandSide : public Function template double PressureRightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { return 0; } @@ -171,11 +171,11 @@ PressureRightHandSide::value (const Point &/*p*/, // mentioned in the introduction, we choose a // linear pressure field: template -class PressureBoundaryValues : public Function +class PressureBoundaryValues : public Function { public: PressureBoundaryValues () : Function(1) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; }; @@ -184,7 +184,7 @@ class PressureBoundaryValues : public Function template double PressureBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { return 1-p[0]; } @@ -201,11 +201,11 @@ PressureBoundaryValues::value (const Point &p, // values. This is as explained in the // introduction: template -class SaturationBoundaryValues : public Function +class SaturationBoundaryValues : public Function { public: SaturationBoundaryValues () : Function(1) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; }; @@ -215,7 +215,7 @@ class SaturationBoundaryValues : public Function template double SaturationBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { if (p[0] == 0) return 1; @@ -245,15 +245,15 @@ SaturationBoundaryValues::value (const Point &p, // later go back and choose a different // function for initial values. template -class InitialValues : public Function +class InitialValues : public Function { public: InitialValues () : Function(dim+2) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; - virtual void vector_value (const Point &p, + virtual void vector_value (const Point &p, Vector &value) const; }; @@ -262,7 +262,7 @@ class InitialValues : public Function template double InitialValues::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { return ZeroFunction(dim+2).value (p, component); } @@ -271,7 +271,7 @@ InitialValues::value (const Point &p, template void InitialValues::vector_value (const Point &p, - Vector &values) const + Vector &values) const { ZeroFunction(dim+2).vector_value (p, values); } @@ -312,7 +312,7 @@ namespace SingleCurvingCrack : TensorFunction<2,dim> () {} - + virtual void value_list (const std::vector > &points, std::vector > &values) const; }; @@ -332,12 +332,12 @@ namespace SingleCurvingCrack const double distance_to_flowline = std::fabs(points[p][1]-0.5-0.1*std::sin(10*points[p][0])); - + const double permeability = std::max(std::exp(-(distance_to_flowline* distance_to_flowline) / (0.1 * 0.1)), 0.01); - + for (unsigned int d=0; d () {} - + virtual void value_list (const std::vector > &points, std::vector > &values) const; @@ -432,7 +432,7 @@ namespace RandomMedium (dim == 3 ? 100 : throw ExcNotImplemented())); - + std::vector > centers_list (N); for (unsigned int i=0; i::InverseMatrix (const Matrix &m) {} - + template void InverseMatrix::vmult (Vector &dst, const Vector &src) const @@ -530,8 +530,8 @@ void InverseMatrix::vmult (Vector &dst, SolverCG<> cg (solver_control); dst = 0; - - cg.solve (*matrix, dst, src, PreconditionIdentity()); + + cg.solve (*matrix, dst, src, PreconditionIdentity()); } @@ -548,7 +548,7 @@ class SchurComplement : public Subscriptor private: const SmartPointer > system_matrix; const SmartPointer > > m_inverse; - + mutable Vector tmp1, tmp2; }; @@ -585,7 +585,7 @@ class ApproximateSchurComplement : public Subscriptor private: const SmartPointer > system_matrix; - + mutable Vector tmp1, tmp2; }; @@ -653,14 +653,14 @@ TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) template void TwoPhaseFlowProblem::make_grid_and_dofs () { - GridGenerator::hyper_cube (triangulation, 0, 1); + GridGenerator::hyper_cube (triangulation, 0, 1); triangulation.refine_global (n_refinement_steps); - - dof_handler.distribute_dofs (fe); + + dof_handler.distribute_dofs (fe); DoFRenumbering::component_wise (dof_handler); - + std::vector dofs_per_component (dim+2); - DoFTools::count_dofs_per_component (dof_handler, dofs_per_component); + DoFTools::count_dofs_per_component (dof_handler, dofs_per_component); const unsigned int n_u = dofs_per_component[0], n_p = dofs_per_component[dim], n_s = dofs_per_component[dim+1]; @@ -673,10 +673,10 @@ void TwoPhaseFlowProblem::make_grid_and_dofs () << " (" << n_u << '+' << n_p << '+'<< n_s <<')' << std::endl << std::endl; - + const unsigned int n_couplings = dof_handler.max_couplings_between_dofs(); - + sparsity_pattern.reinit (3,3); sparsity_pattern.block(0,0).reinit (n_u, n_u, n_couplings); sparsity_pattern.block(1,0).reinit (n_p, n_u, n_couplings); @@ -687,28 +687,28 @@ void TwoPhaseFlowProblem::make_grid_and_dofs () sparsity_pattern.block(0,2).reinit (n_u, n_s, n_couplings); sparsity_pattern.block(1,2).reinit (n_p, n_s, n_couplings); sparsity_pattern.block(2,2).reinit (n_s, n_s, n_couplings); - + sparsity_pattern.collect_sizes(); - + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress(); - + system_matrix.reinit (sparsity_pattern); - + solution.reinit (3); solution.block(0).reinit (n_u); solution.block(1).reinit (n_p); solution.block(2).reinit (n_s); solution.collect_sizes (); - + old_solution.reinit (3); old_solution.block(0).reinit (n_u); old_solution.block(1).reinit (n_p); old_solution.block(2).reinit (n_s); old_solution.collect_sizes (); - + system_rhs.reinit (3); system_rhs.block(0).reinit (n_u); system_rhs.block(1).reinit (n_p); @@ -742,23 +742,23 @@ void TwoPhaseFlowProblem::make_grid_and_dofs () // crack permeability function is as simple // as just changing the namespace name. template -void TwoPhaseFlowProblem::assemble_system () -{ +void TwoPhaseFlowProblem::assemble_system () +{ system_matrix=0; system_rhs=0; - QGauss quadrature_formula(degree+2); + QGauss quadrature_formula(degree+2); QGauss face_quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, + FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); - FEFaceValues fe_face_values (fe, face_quadrature_formula, + FEFaceValues fe_face_values (fe, face_quadrature_formula, update_values | update_normal_vectors | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; - + const unsigned int n_q_points = quadrature_formula.size(); const unsigned int n_face_q_points = face_quadrature_formula.size(); @@ -766,15 +766,15 @@ void TwoPhaseFlowProblem::assemble_system () Vector local_rhs (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - + const PressureRightHandSide pressure_right_hand_side; const PressureBoundaryValues pressure_boundary_values; - const RandomMedium::KInverse k_inverse; - + const RandomMedium::KInverse k_inverse; + std::vector pressure_rhs_values (n_q_points); std::vector boundary_values (n_face_q_points); std::vector > k_inverse_values (n_q_points); - + std::vector > old_solution_values(n_q_points, Vector(dim+2)); std::vector > > old_solution_grads(n_q_points, std::vector > (dim+2)); @@ -787,7 +787,7 @@ void TwoPhaseFlowProblem::assemble_system () cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) - { + { fe_values.reinit (cell); local_matrix = 0; local_rhs = 0; @@ -830,7 +830,7 @@ void TwoPhaseFlowProblem::assemble_system () // self-explanatory given the explicit // form of the bilinear form stated in // the introduction: - for (unsigned int q=0; q::assemble_system () const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q); const double div_phi_i_u = fe_values[velocities].divergence (i, q); const double phi_i_p = fe_values[pressure].value (i, q); - const double phi_i_s = fe_values[saturation].value (i, q); - const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient(i, q); - + const double phi_i_s = fe_values[saturation].value (i, q); + for (unsigned int j=0; j phi_j_u = fe_values[velocities].value (j, q); const double div_phi_j_u = fe_values[velocities].divergence (j, q); const double phi_j_p = fe_values[pressure].value (j, q); const double phi_j_s = fe_values[saturation].value (j, q); - + local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * mobility_inverse(old_s,viscosity) * phi_j_u - div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u + phi_i_s * phi_j_s) - * fe_values.JxW(q); + * fe_values.JxW(q); } local_rhs(i) += (-phi_i_p * pressure_rhs_values[q])* fe_values.JxW(q); } - + // Next, we also have to deal with the // pressure boundary values. This, @@ -870,12 +869,12 @@ void TwoPhaseFlowProblem::assemble_system () if (cell->at_boundary(face_no)) { fe_face_values.reinit (cell, face_no); - + pressure_boundary_values .value_list (fe_face_values.get_quadrature_points(), boundary_values); - for (unsigned int q=0; q @@ -899,7 +898,7 @@ void TwoPhaseFlowProblem::assemble_system () system_matrix.add (local_dof_indices[i], local_dof_indices[j], local_matrix(i,j)); - + for (unsigned int i=0; i::assemble_system () // been computed. We therefore have this // separate function to this end. template -void TwoPhaseFlowProblem::assemble_rhs_S () -{ - QGauss quadrature_formula(degree+2); - QGauss face_quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, +void TwoPhaseFlowProblem::assemble_rhs_S () +{ + QGauss quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); + FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); - FEFaceValues fe_face_values (fe, face_quadrature_formula, + FEFaceValues fe_face_values (fe, face_quadrature_formula, update_values | update_normal_vectors | update_quadrature_points | update_JxW_values); - FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, + FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, update_values); - + const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); const unsigned int n_face_q_points = face_quadrature_formula.size(); - + Vector local_rhs (dofs_per_cell); std::vector > old_solution_values(n_q_points, Vector(dim+2)); @@ -952,7 +951,7 @@ void TwoPhaseFlowProblem::assemble_rhs_S () SaturationBoundaryValues saturation_boundary_values; const FEValuesExtractors::Scalar saturation (dim+1); - + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -970,7 +969,7 @@ void TwoPhaseFlowProblem::assemble_rhs_S () // \mathbf{v}^{n+1},\nabla sigma)$, // where $\sigma$ is the saturation // component of the test function: - for (unsigned int q=0; q::assemble_rhs_S () const double phi_i_s = fe_values[saturation].value (i, q); const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient (i, q); - + local_rhs(i) += (time_step * f_saturation(old_s,viscosity) * present_u * @@ -1027,15 +1026,15 @@ void TwoPhaseFlowProblem::assemble_rhs_S () neighbor_face = cell->neighbor_of_neighbor(face_no); fe_face_values_neighbor.reinit (neighbor, neighbor_face); - + fe_face_values_neighbor .get_function_values (old_solution, old_solution_values_face_neighbor); - + for (unsigned int q=0; q::assemble_rhs_S () fe_face_values.normal_vector(q); const bool is_outflow_q_point = (normal_flux >= 0); - + for (unsigned int i=0; i::assemble_rhs_S () fe_face_values.JxW(q); } } - + cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i::assemble_rhs_S () // step-20. After that, we have to deal with // the saturation equation (see below): template -void TwoPhaseFlowProblem::solve () +void TwoPhaseFlowProblem::solve () { const InverseMatrix > m_inverse (system_matrix.block(0,0)); Vector tmp (solution.block(0).size()); Vector schur_rhs (solution.block(1).size()); Vector tmp2 (solution.block(2).size()); - + // First the pressure, using the pressure // Schur complement of the first two @@ -1095,24 +1094,24 @@ void TwoPhaseFlowProblem::solve () system_matrix.block(1,0).vmult (schur_rhs, tmp); schur_rhs -= system_rhs.block(1); - + SchurComplement schur_complement (system_matrix, m_inverse); - + ApproximateSchurComplement approximate_schur_complement (system_matrix); - + InverseMatrix preconditioner (approximate_schur_complement); - + SolverControl solver_control (solution.block(1).size(), 1e-12*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); cg.solve (schur_complement, solution.block(1), schur_rhs, preconditioner); - + std::cout << " " << solver_control.last_step() << " CG Schur complement iterations for pressure." @@ -1124,7 +1123,7 @@ void TwoPhaseFlowProblem::solve () system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; tmp += system_rhs.block(0); - + m_inverse.vmult (solution.block(0), tmp); } @@ -1160,25 +1159,25 @@ void TwoPhaseFlowProblem::solve () // reasonable range: assemble_rhs_S (); { - + SolverControl solver_control (system_matrix.block(2,2).m(), 1e-8*system_rhs.block(2).l2_norm()); SolverCG<> cg (solver_control); cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), PreconditionIdentity()); - + project_back_saturation (); - + std::cout << " " << solver_control.last_step() << " CG iterations for saturation." - << std::endl; - } + << std::endl; + } + - - old_solution = solution; + old_solution = solution; } - + // @sect4{TwoPhaseFlowProblem::output_results} @@ -1191,7 +1190,7 @@ void TwoPhaseFlowProblem::output_results () const { if (timestep_number % 5 != 0) return; - + std::vector solution_names; switch (dim) { @@ -1201,7 +1200,7 @@ void TwoPhaseFlowProblem::output_results () const solution_names.push_back ("p"); solution_names.push_back ("S"); break; - + case 3: solution_names.push_back ("u"); solution_names.push_back ("v"); @@ -1209,18 +1208,18 @@ void TwoPhaseFlowProblem::output_results () const solution_names.push_back ("p"); solution_names.push_back ("S"); break; - + default: Assert (false, ExcNotImplemented()); } - + DataOut data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, solution_names); data_out.build_patches (degree+1); - + std::ostringstream filename; filename << "solution-" << timestep_number << ".vtk"; @@ -1281,16 +1280,16 @@ template double TwoPhaseFlowProblem::get_maximal_velocity () const { - QGauss quadrature_formula(degree+2); + QGauss quadrature_formula(degree+2); const unsigned int n_q_points = quadrature_formula.size(); - FEValues fe_values (fe, quadrature_formula, + FEValues fe_values (fe, quadrature_formula, update_values); std::vector > solution_values(n_q_points, Vector(dim+2)); double max_velocity = 0; - + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -1303,8 +1302,8 @@ TwoPhaseFlowProblem::get_maximal_velocity () const { Tensor<1,dim> velocity; for (unsigned int i=0; i::get_maximal_velocity () const // end time of a time step only at the end of // the time step. template -void TwoPhaseFlowProblem::run () +void TwoPhaseFlowProblem::run () { make_grid_and_dofs(); - + { ConstraintMatrix constraints; constraints.close(); - + VectorTools::project (dof_handler, constraints, QGauss(degree+2), InitialValues(), old_solution); } - + timestep_number = 1; double time = 0; - + do - { + { std::cout << "Timestep " << timestep_number - << std::endl; + << std::endl; assemble_system (); solve (); - + output_results (); time += time_step; @@ -1386,7 +1385,7 @@ void TwoPhaseFlowProblem::run () while (time <= 250); } - + // @sect3{The main function} // That's it. In the main function, we pass @@ -1396,7 +1395,7 @@ void TwoPhaseFlowProblem::run () // elements, i.e. $RT_0\times DQ_0 \times // DQ_0$. The rest is as in all the other // programs. -int main () +int main () { try { @@ -1415,10 +1414,10 @@ int main () << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; - + return 1; } - catch (...) + catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" -- 2.39.5