From: Martin Kronbichler Date: Wed, 16 Apr 2008 06:43:10 +0000 (+0000) Subject: Block Schur preconditioned solver on u and p using the BlockMatrixArray class. X-Git-Tag: v8.0.0~9215 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=95251ebb0820811c46a5ed0338287d8825217e1c;p=dealii.git Block Schur preconditioned solver on u and p using the BlockMatrixArray class. git-svn-id: https://svn.dealii.org/trunk@15997 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 5bf011105f..033d8a4830 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -21,9 +21,12 @@ #include #include #include +#include #include #include #include +#include +#include #include #include @@ -57,14 +60,29 @@ using namespace dealii; - template -class BoussinesqFlowProblem + struct InnerPreconditioner; + + template <> + struct InnerPreconditioner<2> + { + typedef SparseDirectUMFPACK type; + }; + + template <> + struct InnerPreconditioner<3> + { + typedef SparseILU type; + }; + + +template +class BoussinesqFlowProblem { public: BoussinesqFlowProblem (const unsigned int degree); void run (); - + private: void setup_dofs (const bool setup_matrices); void assemble_system (); @@ -73,26 +91,27 @@ class BoussinesqFlowProblem void solve (); void output_results () const; void refine_mesh (); - + const unsigned int degree; - + Triangulation triangulation; FESystem fe; DoFHandler dof_handler; ConstraintMatrix hanging_node_constraints; - + BlockSparsityPattern sparsity_pattern; BlockSparseMatrix system_matrix; double time_step; unsigned int timestep_number; - + BlockVector solution; BlockVector old_solution; BlockVector system_rhs; - boost::shared_ptr A_preconditioner; + boost::shared_ptr::type> A_preconditioner; + boost::shared_ptr > Mp_preconditioner; bool rebuild_matrices; bool rebuild_preconditioner; @@ -103,11 +122,11 @@ class BoussinesqFlowProblem 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; }; @@ -116,7 +135,7 @@ class PressureBoundaryValues : public Function template double PressureBoundaryValues::value (const Point &/*p*/, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { return 0; } @@ -124,11 +143,11 @@ PressureBoundaryValues::value (const Point &/*p*/, template -class TemperatureBoundaryValues : public Function +class TemperatureBoundaryValues : public Function { public: TemperatureBoundaryValues () : Function(1) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; }; @@ -138,7 +157,7 @@ class TemperatureBoundaryValues : public Function template double TemperatureBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { //TODO: leftover from olden times. replace by something sensible once we have //diffusion in the temperature field @@ -152,15 +171,15 @@ TemperatureBoundaryValues::value (const Point &p, 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; }; @@ -168,7 +187,7 @@ class InitialValues : public Function template double InitialValues::value (const Point &, - const unsigned int) const + const unsigned int) const { return 0; } @@ -177,7 +196,7 @@ InitialValues::value (const Point &, template void InitialValues::vector_value (const Point &p, - Vector &values) const + Vector &values) const { for (unsigned int c=0; cn_components; ++c) values(c) = InitialValues::value (p, c); @@ -186,15 +205,15 @@ InitialValues::vector_value (const Point &p, template -class RightHandSide : public Function +class RightHandSide : public Function { public: RightHandSide () : 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; }; @@ -202,7 +221,7 @@ class RightHandSide : public Function template double RightHandSide::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { if (component == dim+1) return ((p.distance (Point(.3,.1)) < 1./32) @@ -222,7 +241,7 @@ RightHandSide::value (const Point &p, template void RightHandSide::vector_value (const Point &p, - Vector &values) const + Vector &values) const { for (unsigned int c=0; cn_components; ++c) values(c) = RightHandSide::value (p, c); @@ -256,7 +275,7 @@ InverseMatrix::InverseMatrix (const Matrix &m, {} - + template void InverseMatrix::vmult (Vector &dst, const Vector &src) const @@ -277,45 +296,67 @@ void InverseMatrix::vmult (Vector &dst, } - -template -class SchurComplement : public Subscriptor + // This is the implementation + // of the Schur complement + // preconditioner as described + // in the section on improved + // solvers in step-22. See there + // for more explanation of the + // method. +template +class BlockSchurPreconditioner : public Subscriptor { public: - SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix,Preconditioner> &Minv); + BlockSchurPreconditioner (const BlockSparseMatrix &S, + const InverseMatrix,PreconditionerMp> &Mpinv, + const PreconditionerA &Apreconditioner); - void vmult (Vector &dst, - const Vector &src) const; + void vmult (BlockVector &dst, + const BlockVector &src) const; private: const SmartPointer > system_matrix; - const SmartPointer,Preconditioner> > m_inverse; - - mutable Vector tmp1, tmp2; -}; + const SmartPointer, + PreconditionerMp > > m_inverse; + const PreconditionerA &a_preconditioner; + mutable Vector tmp; +}; -template -SchurComplement:: -SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix,Preconditioner> &Minv) +template +BlockSchurPreconditioner::BlockSchurPreconditioner( + const BlockSparseMatrix &S, + const InverseMatrix,PreconditionerMp> &Mpinv, + const PreconditionerA &Apreconditioner + ) : - system_matrix (&A), - m_inverse (&Minv), - tmp1 (A.block(0,0).m()), - tmp2 (A.block(0,0).m()) -{} - + system_matrix (&S), + m_inverse (&Mpinv), + a_preconditioner (Apreconditioner), + tmp (S.block(1,1).m()) +{ +} -template -void SchurComplement::vmult (Vector &dst, - const Vector &src) const + // Now the interesting function, the multiplication of + // the preconditioner with a BlockVector. +template +void BlockSchurPreconditioner::vmult ( + BlockVector &dst, + const BlockVector &src) const { - system_matrix->block(0,1).vmult (tmp1, src); - m_inverse->vmult (tmp2, tmp1); - system_matrix->block(1,0).vmult (dst, tmp2); + // Form u_new = A^{-1} u + a_preconditioner.vmult (dst.block(0), src.block(0)); + // Form tmp = - B u_new + p + // (SparseMatrix::residual + // does precisely this) + system_matrix->block(1,0).residual(tmp, + dst.block(0), src.block(1)); + // Change sign in tmp.block(1) + tmp *= -1; + // Multiply by approximate Schur complement + // (i.e. a pressure mass matrix) + m_inverse->vmult (dst.block(1), tmp); } @@ -338,8 +379,8 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) template void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) -{ - dof_handler.distribute_dofs (fe); +{ + dof_handler.distribute_dofs (fe); DoFRenumbering::component_wise (dof_handler); hanging_node_constraints.clear (); @@ -353,7 +394,7 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) hanging_node_constraints.close (); 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] * dim, n_p = dofs_per_component[dim], n_T = dofs_per_component[dim+1]; @@ -366,14 +407,14 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) << " (" << n_u << '+' << n_p << '+'<< n_T <<')' << std::endl << std::endl; - + const unsigned int n_couplings = dof_handler.max_couplings_between_dofs(); if (setup_matrices == true) { system_matrix.clear (); - + 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); @@ -384,29 +425,29 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) sparsity_pattern.block(0,2).reinit (n_u, n_T, n_couplings); sparsity_pattern.block(1,2).reinit (n_p, n_T, n_couplings); sparsity_pattern.block(2,2).reinit (n_T, n_T, n_couplings); - - sparsity_pattern.collect_sizes(); - + sparsity_pattern.collect_sizes(); + + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); hanging_node_constraints.condense (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_T); 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_T); old_solution.collect_sizes (); - + system_rhs.reinit (3); system_rhs.block(0).reinit (n_u); system_rhs.block(1).reinit (n_p); @@ -417,14 +458,14 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) template -void BoussinesqFlowProblem::assemble_system () +void BoussinesqFlowProblem::assemble_system () { if (rebuild_matrices == true) 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, @@ -436,12 +477,12 @@ void BoussinesqFlowProblem::assemble_system () update_gradients : UpdateFlags(0))); - 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(); @@ -449,17 +490,17 @@ void BoussinesqFlowProblem::assemble_system () Vector local_rhs (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - + const PressureBoundaryValues pressure_boundary_values; - + std::vector boundary_values (n_face_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)); const double Raleigh_number = 10; - + const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); const FEValuesExtractors::Scalar temperature (dim+1); @@ -468,7 +509,7 @@ void BoussinesqFlowProblem::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; @@ -478,7 +519,7 @@ void BoussinesqFlowProblem::assemble_system () 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_T = fe_values[temperature].value (i, q); + const double phi_i_T = fe_values[temperature].value (i, q); const Tensor<1,dim> grad_phi_i_T = fe_values[temperature].gradient(i, q); - + for (unsigned int j=0; j @@ -501,24 +542,24 @@ void BoussinesqFlowProblem::assemble_system () 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_T = fe_values[temperature].value (j, q); - + local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u - div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u + phi_i_p * phi_j_p + phi_i_T * phi_j_T) - * fe_values.JxW(q); + * fe_values.JxW(q); } } - + const Point gravity (0,1); - + local_rhs(i) += (Raleigh_number * gravity * phi_i_u * old_temperature)* fe_values.JxW(q); } } - + for (unsigned int face_no=0; face_no::faces_per_cell; @@ -526,12 +567,12 @@ void BoussinesqFlowProblem::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 @@ -554,16 +595,16 @@ void BoussinesqFlowProblem::assemble_system () local_dof_indices[j], local_matrix(i,j)); } - + for (unsigned int i=0; i boundary_values; @@ -581,11 +622,11 @@ void BoussinesqFlowProblem::assemble_system () // { // std::cout << "Found cell and vertex: " << cell << ' ' // << v << std::endl; - + // boundary_values[cell->vertex_dof_index(v,0)] = 0; // break; -// } - +// } + // std::vector component_mask (dim+2, true); // component_mask[dim] = component_mask[dim+1] = false; // VectorTools::interpolate_boundary_values (dof_handler, @@ -597,7 +638,7 @@ void BoussinesqFlowProblem::assemble_system () // MatrixTools::apply_boundary_values (boundary_values, // system_matrix, // solution, -// system_rhs); +// system_rhs); } if (rebuild_preconditioner == true) @@ -605,15 +646,29 @@ void BoussinesqFlowProblem::assemble_system () Assert (rebuild_matrices == true, ExcMessage ("There is no point in rebuilding the preconditioner " "without a rebuilt matrix!")); - + std::cout << " Rebuilding preconditioner..." << std::flush; - + + // Rebuild the preconditioner + // for the velocity-velocity + // block (0,0) A_preconditioner - = boost::shared_ptr(new SparseDirectUMFPACK()); - A_preconditioner->initialize (system_matrix.block(0,0)); + = boost::shared_ptr::type> + (new typename InnerPreconditioner::type()); + A_preconditioner->initialize (system_matrix.block(0,0), + typename InnerPreconditioner::type::AdditionalData()); + + // Rebuild the preconditioner + // for the pressure-pressure + // block (1,1) + Mp_preconditioner + = boost::shared_ptr > + (new SparseILU); + Mp_preconditioner->initialize (system_matrix.block(0,0), + SparseILU::AdditionalData()); std::cout << std::endl; - + rebuild_preconditioner = false; } @@ -626,28 +681,28 @@ void BoussinesqFlowProblem::assemble_system () template -void BoussinesqFlowProblem::assemble_rhs_T () -{ - QGauss quadrature_formula(degree+2); +void BoussinesqFlowProblem::assemble_rhs_T () +{ + 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); - FESubfaceValues fe_subface_values (fe, face_quadrature_formula, + FESubfaceValues fe_subface_values (fe, face_quadrature_formula, update_values | update_normal_vectors | update_JxW_values); - FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, + FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, update_values); - FESubfaceValues fe_subface_values_neighbor (fe, face_quadrature_formula, + FESubfaceValues fe_subface_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)); @@ -660,14 +715,14 @@ void BoussinesqFlowProblem::assemble_rhs_T () std::vector > > present_solution_grads(n_q_points, std::vector >(dim+2)); - - + + std::vector neighbor_temperature (n_face_q_points); std::vector local_dof_indices (dofs_per_cell); TemperatureBoundaryValues temperature_boundary_values; const FEValuesExtractors::Scalar temperature (dim+1); - + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -681,7 +736,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_values.get_function_values (solution, present_solution_values); fe_values.get_function_gradients (solution, present_solution_grads); - for (unsigned int q=0; q::assemble_rhs_T () double present_div_u = 0; for (unsigned int d=0; d grad_phi_i_T = fe_values[temperature].gradient (i, q); const Point p = fe_values.quadrature_point(q); - + local_rhs(i) += (time_step * old_T * (present_u * @@ -715,8 +770,8 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_values.JxW(q); } - -//TODO: unify the code that actually does the assembly down below + +//TODO: unify the code that actually does the assembly down below for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) if (cell->at_boundary(face_no) @@ -746,14 +801,14 @@ void BoussinesqFlowProblem::assemble_rhs_T () { const typename DoFHandler::active_cell_iterator neighbor = cell->neighbor(face_no); - + fe_face_values_neighbor.reinit (neighbor, cell->neighbor_of_neighbor(face_no)); - + fe_face_values_neighbor .get_function_values (old_solution, old_solution_values_face_neighbor); - + for (unsigned int q=0; q::assemble_rhs_T () fe_face_values.normal_vector(q); const bool is_outflow_q_point = (normal_flux >= 0); - + for (unsigned int i=0; i::assemble_rhs_T () fe_face_values_neighbor.reinit (neighbor, cell->neighbor_of_neighbor(face_no)); - + fe_face_values_neighbor .get_function_values (old_solution, old_solution_values_face_neighbor); - + for (unsigned int q=0; q::assemble_rhs_T () fe_subface_values.normal_vector(q); const bool is_outflow_q_point = (normal_flux >= 0); - + for (unsigned int i=0; i::assemble_rhs_T () fe_subface_values_neighbor.reinit (neighbor, neighbor_face_no, neighbor_subface_no); - + fe_subface_values_neighbor .get_function_values (old_solution, old_solution_values_face_neighbor); - + for (unsigned int q=0; q::assemble_rhs_T () fe_face_values.normal_vector(q); const bool is_outflow_q_point = (normal_flux >= 0); - + for (unsigned int i=0; i::assemble_rhs_T () fe_face_values[temperature].value (i,q) * fe_face_values.JxW(q); } - } - + } + cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i -void BoussinesqFlowProblem::solve () +void BoussinesqFlowProblem::solve () { solution = old_solution; - - const InverseMatrix,SparseDirectUMFPACK> - A_inverse (system_matrix.block(0,0), *A_preconditioner); - Vector tmp (solution.block(0).size()); - Vector schur_rhs (solution.block(1).size()); - Vector tmp2 (solution.block(2).size()); - + + // Use the BlockMatrixArray structure + // for extracting only the upper left + // 2x2 blocks from the matrix that will + // be used for the solution of the + // blocked system. { - A_inverse.vmult (tmp, system_rhs.block(0)); - system_matrix.block(1,0).vmult (schur_rhs, tmp); - schur_rhs -= system_rhs.block(1); - - - SchurComplement - schur_complement (system_matrix, A_inverse); - - SolverControl solver_control (system_matrix.block(0,0).m(), - 1e-6*schur_rhs.l2_norm()); - SolverCG<> cg (solver_control); - - PreconditionSSOR<> preconditioner; - preconditioner.initialize (system_matrix.block(1,1), 1.2); - - InverseMatrix,PreconditionSSOR<> > - m_inverse (system_matrix.block(1,1), preconditioner); - - try - { - cg.solve (schur_complement, solution.block(1), schur_rhs, - m_inverse); - } - catch (...) - { - abort (); - } + GrowingVectorMemory > simple_mem; + BlockMatrixArray stokes_submatrix(2, 2, simple_mem); - // produce a consistent flow field - hanging_node_constraints.distribute (solution); - - std::cout << " " - << solver_control.last_step() - << " CG Schur complement iterations for pressure." - << std::endl; - } + stokes_submatrix.enter(system_matrix.block(0,0),0,0); + stokes_submatrix.enter(system_matrix.block(0,1),0,1); + stokes_submatrix.enter(system_matrix.block(1,0),1,0); - { - system_matrix.block(0,1).vmult (tmp, solution.block(1)); - tmp *= -1; - tmp += system_rhs.block(0); - - A_inverse.vmult (solution.block(0), tmp); + // Define some temporary vectors + // for the solution process. + std::vector block_sizes(2); + block_sizes[0] = solution.block(0).size(); + block_sizes[1] = solution.block(1).size(); - // produce a consistent pressure field - hanging_node_constraints.distribute (solution); - } + BlockVector up_rhs(block_sizes); + BlockVector up(block_sizes); + + up_rhs.block(0) = system_rhs.block(0); + up_rhs.block(1) = system_rhs.block(1); + + // Set up inverse matrix for + // pressure mass matrix + InverseMatrix,SparseILU > + mp_inverse (system_matrix.block(1,1), *Mp_preconditioner); + + // Set up block Schur preconditioner + BlockSchurPreconditioner::type, + SparseILU > + preconditioner (system_matrix, mp_inverse, *A_preconditioner); + + // Set up GMRES solver and + // solve. + SolverControl solver_control (system_matrix.m(), + 1e-6*system_rhs.l2_norm()); + SolverGMRES > gmres(solver_control, + SolverGMRES >::AdditionalData(100)); + + gmres.solve(stokes_submatrix, up, up_rhs, + preconditioner); + + // Produce a constistent solution field + hanging_node_constraints.distribute (up); + + solution.block(0) = up.block(0); + solution.block(1) = up.block(1); + } // for DGQ1 needs to be /15 time_step = GridTools::minimal_cell_diameter(triangulation) / std::max (get_maximal_velocity(), .05) / 2; assemble_rhs_T (); { - + SolverControl solver_control (system_matrix.block(2,2).m(), 1e-8*system_rhs.block(2).l2_norm()); SolverCG<> cg (solver_control); PreconditionJacobi<> preconditioner; preconditioner.initialize (system_matrix.block(2,2)); - + try { cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), @@ -986,10 +1038,10 @@ void BoussinesqFlowProblem::solve () { abort (); } - + // produce a consistent temperature field hanging_node_constraints.distribute (solution); - + std::cout << " " << solver_control.last_step() << " CG iterations for temperature." @@ -998,9 +1050,9 @@ void BoussinesqFlowProblem::solve () << *std::max_element (solution.block(2).begin(), solution.block(2).end()) << std::endl; - } + } } - + template @@ -1008,11 +1060,11 @@ void BoussinesqFlowProblem::output_results () const { if (timestep_number % 10 != 0) return; - + std::vector solution_names (dim, "velocity"); solution_names.push_back ("p"); solution_names.push_back ("T"); - + DataOut data_out; data_out.attach_dof_handler (dof_handler); @@ -1023,13 +1075,13 @@ void BoussinesqFlowProblem::output_results () const for (unsigned int i=0; i::type_dof_data, data_component_interpretation); - + data_out.build_patches (); - + std::ostringstream filename; filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk"; @@ -1041,11 +1093,11 @@ void BoussinesqFlowProblem::output_results () const template void -BoussinesqFlowProblem::refine_mesh () +BoussinesqFlowProblem::refine_mesh () { Vector estimated_error_per_cell (triangulation.n_active_cells()); -//TODO do this better +//TODO do this better DerivativeApproximation::approximate_gradient (dof_handler, old_solution, estimated_error_per_cell, @@ -1069,7 +1121,7 @@ BoussinesqFlowProblem::refine_mesh () Vector x_old_solution (dof_handler.n_dofs()); x_old_solution = old_solution; - + soltrans.prepare_for_coarsening_and_refinement(x_old_solution); triangulation.execute_coarsening_and_refinement (); @@ -1090,16 +1142,16 @@ template double BoussinesqFlowProblem::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(); @@ -1112,20 +1164,20 @@ BoussinesqFlowProblem::get_maximal_velocity () const { Tensor<1,dim> velocity; for (unsigned int i=0; i -void BoussinesqFlowProblem::run () +void BoussinesqFlowProblem::run () { switch (dim) { @@ -1137,7 +1189,7 @@ void BoussinesqFlowProblem::run () // triangulation.set_boundary (0, boundary); GridGenerator::hyper_cube (triangulation); - + triangulation.refine_global (6); break; @@ -1147,10 +1199,10 @@ void BoussinesqFlowProblem::run () { GridGenerator::hyper_shell (triangulation, Point(), 0.5, 1.0); - + static HyperShellBoundary boundary; triangulation.set_boundary (0, boundary); - + triangulation.refine_global (2); break; @@ -1159,7 +1211,7 @@ void BoussinesqFlowProblem::run () default: Assert (false, ExcNotImplemented()); } - + const bool do_adaptivity = false; @@ -1172,7 +1224,7 @@ void BoussinesqFlowProblem::run () QGauss(degree+2), InitialValues(), old_solution); - + for (unsigned int pre_refinement=0; pre_refinement<4-dim; ++pre_refinement) { refine_mesh (); @@ -1194,29 +1246,29 @@ void BoussinesqFlowProblem::run () InitialValues(), old_solution); } - + timestep_number = 0; double time = 0; - + do - { + { std::cout << "Timestep " << timestep_number << ": t=" << time << ", dt=" << time_step - << std::endl; + << std::endl; std::cout << " Assembling..." << std::endl; - assemble_system (); + assemble_system (); std::cout << " Solving..." << std::endl; solve (); - + output_results (); time += time_step; ++timestep_number; - - old_solution = solution; + + old_solution = solution; std::cout << std::endl; @@ -1227,9 +1279,9 @@ void BoussinesqFlowProblem::run () while (time <= 5); } - -int main () + +int main () { try { @@ -1248,10 +1300,10 @@ int main () << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; - + return 1; } - catch (...) + catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------"