From 939f1320997d0c805ff8001a815c9e5bc6fca838 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 23 Feb 2008 19:26:27 +0000 Subject: [PATCH] Re-indent the entire file. Fix up a few places with formulas and other formatting. git-svn-id: https://svn.dealii.org/trunk@15760 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 1066 +++++++++++++-------------- 1 file changed, 517 insertions(+), 549 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index a726cb4ed3..516dfbb38b 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -12,10 +12,10 @@ /* further information on this license. */ - // @sect3{Include files} + // @sect3{Include files} - // As usual, we start by including - // some well-known files. + // As usual, we start by including + // some well-known files. #include #include #include @@ -51,53 +51,51 @@ #include #include - // As in - // step-29, we include the libary - // for the sparse direct solver - // UMFPACK. + // As in step-29, we include the libary for + // the sparse direct solver UMFPACK. #include - // This includes the libary for the - // incomplete LU factorization that will - // be used as a preconditioner in 3D. + // This includes the libary for the + // incomplete LU factorization that will + // be used as a preconditioner in 3D. #include - // This is C++: + // This is C++: #include #include - // As in all programs, the namespace is set - // to dealii. + // As in all programs, the namespace is set + // to dealii. using namespace dealii; - // @sect3{Defining the inner preconditioner type} + // @sect3{Defining the inner preconditioner type} - // As explained in the introduction, we - // are going to use different preconditioners - // for two and three space dimensions, - // respectively. We differentiate between - // them by the use of the spatial dimension - // as a template parameter. See step-4 for - // details on templates. - // We are not going to create any preconditioner - // object here, all we do is to create a - // data structure that holds the information - // on it so we can write our program in a - // dimension-independent way. + // As explained in the introduction, we are + // going to use different preconditioners for + // two and three space dimensions, + // respectively. We differentiate between + // them by the use of the spatial dimension + // as a template parameter. See step-4 for + // details on templates. We are not going to + // create any preconditioner object here, all + // we do is to create a data structure that + // holds the information on it so we can + // write our program in a + // dimension-independent way. template struct InnerPreconditioner; - // In 2D, we are going to use a sparse direct - // solve as preconditioner. The syntax is - // known from step-29. + // In 2D, we are going to use a sparse direct + // solve as preconditioner. The syntax is + // known from step-29. template <> struct InnerPreconditioner<2> { typedef SparseDirectUMFPACK type; }; - // And the ILU preconditioning in 3D, called - // by SparseILU@. + // And the ILU preconditioning in 3D, called + // by SparseILU@. template <> struct InnerPreconditioner<3> { @@ -105,14 +103,13 @@ struct InnerPreconditioner<3> }; - // @sect3{The StokesProblem class template} + // @sect3{The StokesProblem class template} - // This is an adaptation of step-20, - // so the main class and the data types - // are the same as used there. In this - // example we also use adaptive grid - // refinement, which is handled in complete - // analogy to step-6. + // This is an adaptation of step-20, so the + // main class and the data types are the same + // as used there. In this example we also use + // adaptive grid refinement, which is handled + // in complete analogy to step-6. template class StokesProblem { @@ -141,38 +138,39 @@ class StokesProblem BlockVector solution; BlockVector system_rhs; - // This one is new: We shall use a so-called - // shared pointer structure to access - // the preconditioner. This provides - // flexibility when using the object - // that the pointer refers to, as e.g. - // the reset option. + // This one is new: We shall use a + // so-called shared pointer structure to + // access the preconditioner. This + // provides flexibility when using the + // object that the pointer refers to, as + // e.g. the reset option. boost::shared_ptr::type> A_preconditioner; }; - // @sect3{Boundary values and right hand side} - - // As in step-20 and most other example - // programs, the next task is to define - // the parameter functions for the PDE: - // For the Stokes problem, we are going to - // use pressure boundary values at some portion - // of the boundary (Neumann-type), and - // boundary conditions on the velocity - // (Dirichlet type) on the rest of the boundary. - // The pressure boundary condition is - // scalar, and so is the respective function, - // whereas the Dirichlet (velocity) - // condition is vector-valued. Due to the - // structure of deal.II's libraries, we have to - // define the function on the (u,p)-space, but - // we are going to filter out the pressure - // component when condensating the Dirichlet - // data in assemble_system. + // @sect3{Boundary values and right hand side} + + // As in step-20 and most other example + // programs, the next task is to define the + // parameter functions for the PDE: For the + // Stokes problem, we are going to use + // pressure boundary values at some portion + // of the boundary (Neumann-type), and + // boundary conditions on the velocity + // (Dirichlet type) on the rest of the + // boundary. The pressure boundary condition + // is scalar, and so is the respective + // function, whereas the Dirichlet (velocity) + // condition is vector-valued. Due to the + // structure of deal.II's libraries, we have + // to define the function on the (u,p)-space, + // but we are going to filter out the + // pressure component when condensating the + // Dirichlet data in + // assemble_system. - // Given the problem described in the - // introduction, we know which values to - // set for the respective functions. + // Given the problem described in the + // introduction, we know which values to + // set for the respective functions. template class PressureBoundaryValues : public Function { @@ -211,7 +209,7 @@ class BoundaryValues : public Function template double BoundaryValues::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { if (component == 0) return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0)); @@ -222,7 +220,7 @@ BoundaryValues::value (const Point &p, template void BoundaryValues::vector_value (const Point &p, - Vector &values) const + Vector &values) const { for (unsigned int c=0; cn_components; ++c) values(c) = BoundaryValues::value (p, c); @@ -230,8 +228,8 @@ BoundaryValues::vector_value (const Point &p, - // We implement similar functions - // for the right hand side. + // We implement similar functions + // for the right hand side. template class RightHandSide : public Function { @@ -266,27 +264,29 @@ RightHandSide::vector_value (const Point &p, } - // @sect3{Linear solvers and preconditioners} + // @sect3{Linear solvers and preconditioners} - // The linear solvers and preconditioners are - // discussed extensively in the introduction. Here, - // we create the respective objects that will be used. + // The linear solvers and preconditioners are + // discussed extensively in the + // introduction. Here, we create the + // respective objects that will be used. - // @sect4{The InverseMatrix class template} + // @sect4{The InverseMatrix class template} - // This is going to represent the data - // structure for an inverse matrix. This class - // is derived from the one in step-20. The - // only difference is that we now - // do include a preconditioner to the matrix. - // This is going to happen via a template parameter - // class Preconditioner, so - // the preconditioner type will be set when - // an InverseMatrix object is - // created. The member function - // vmult is, as in - // step-20, a multiplication with a vector, - // obtained by solving a linear system. + // This is going to represent the data + // structure for an inverse matrix. This + // class is derived from the one in + // step-20. The only difference is that we + // now do include a preconditioner to the + // matrix. This is going to happen via a + // template parameter class + // Preconditioner, so the + // preconditioner type will be set when an + // InverseMatrix object is + // created. The member function + // vmult is, as in step-20, a + // multiplication with a vector, obtained by + // solving a linear system. template class InverseMatrix : public Subscriptor { @@ -307,62 +307,57 @@ class InverseMatrix : public Subscriptor template InverseMatrix::InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner) - : - matrix (&m), - preconditioner (preconditioner) + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (preconditioner) {} - // This is the implementation of the - // vmult function. We note - // two things: + // This is the implementation of the + // vmult function. We note + // two things: - // Firstly, we use - // a rather large tolerance for the - // solver control. The reason for this - // is that the function is used very - // frequently, and hence, any additional - // effort to make the residual in - // the CG solve smaller makes the - // solution more expensive. Note that - // we do not only use this class as a - // preconditioner for the Schur complement, - // but also when forming the inverse of - // the Laplace matrix - which has to - // be accurate in order to obtain a - // solution to the right problem. + // Firstly, we use a rather large tolerance + // for the solver control. The reason for + // this is that the function is used very + // frequently, and hence, any additional + // effort to make the residual in the CG + // solve smaller makes the solution more + // expensive. Note that we do not only use + // this class as a preconditioner for the + // Schur complement, but also when forming + // the inverse of the Laplace matrix - which + // has to be accurate in order to obtain a + // solution to the right problem. - // Secondly, we catch exceptions from - // the solver at this stage. While this - // is not of greater interest our - // general setting with the requirement - // of accurate inverses (and we indeed - // abort the program when any exception - // occurs), the situation would - // change if an object of the class - // InverseMatrix is only - // used for preconditioning. In such a - // setting, one could imagine to use - // a few CG sweeps as a preconditioner - - // which is done e.g. for mass - // matrices, see the results section - // below. Using - // catch (SolverControl::NoConvergence) {} - // in conjunction with only a few iterations, - // say 10, would result in that effect - - // the program would continue to run - // even though the solver has not converged. - // Note, though, that applying the CG method - // is not a linear operation (see the - // actual CG algorithm for details - // on that), so unconverged - // preconditioners are to be used with - // care in order to not yield a wrong - // solution. + // Secondly, we catch exceptions from the + // solver at this stage. While this is not of + // greater interest our general setting with + // the requirement of accurate inverses (and + // we indeed abort the program when any + // exception occurs), the situation would + // change if an object of the class + // InverseMatrix is only used + // for preconditioning. In such a setting, + // one could imagine to use a few CG sweeps + // as a preconditioner - which is done + // e.g. for mass matrices, see the results + // section below. Using catch + // (SolverControl::NoConvergence) {} + // in conjunction with only a few iterations, + // say 10, would result in that effect - the + // program would continue to run even though + // the solver has not converged. Note, + // though, that applying the CG method is not + // a linear operation (see the actual CG + // algorithm for details on that), so + // unconverged preconditioners are to be used + // with care in order to not yield a wrong + // solution. template void InverseMatrix::vmult (Vector &dst, - const Vector &src) const + const Vector &src) const { SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); SolverCG<> cg (solver_control, vector_memory); @@ -380,31 +375,30 @@ void InverseMatrix::vmult (Vector &dst, } - // @sect4{The SchurComplement class template} - - // This class implements the Schur - // complement discussed in the introduction. - // It is in analogy to step-20. - // Though, we now call it with a template - // parameter Preconditioner - // in order to access that when specifying - // the respective type of the inverse - // matrix class. As a consequence of the - // definition above, the declaration - // InverseMatrix - // now contains the second template - // parameter from preconditioning as above, - // which effects the SmartPointer - // object m_inverse as well. + // @sect4{The SchurComplement class template} + + // This class implements the Schur complement + // discussed in the introduction. It is in + // analogy to step-20. Though, we now call + // it with a template parameter + // Preconditioner in order to + // access that when specifying the respective + // type of the inverse matrix class. As a + // consequence of the definition above, the + // declaration InverseMatrix now + // contains the second template parameter + // from preconditioning as above, which + // effects the SmartPointer + // object m_inverse as well. template class SchurComplement : public Subscriptor { public: SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix, Preconditioner> &Minv); + const InverseMatrix, Preconditioner> &Minv); void vmult (Vector &dst, - const Vector &src) const; + const Vector &src) const; private: const SmartPointer > system_matrix; @@ -418,18 +412,18 @@ class SchurComplement : public Subscriptor template SchurComplement:: SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix,Preconditioner> &Minv) - : - system_matrix (&A), - m_inverse (&Minv), - tmp1 (A.block(0,0).m()), - tmp2 (A.block(0,0).m()) + const InverseMatrix,Preconditioner> &Minv) + : + system_matrix (&A), + m_inverse (&Minv), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) {} template void SchurComplement::vmult (Vector &dst, - const Vector &src) const + const Vector &src) const { system_matrix->block(0,1).vmult (tmp1, src); m_inverse->vmult (tmp2, tmp1); @@ -437,30 +431,32 @@ void SchurComplement::vmult (Vector &dst, } - // @sect3{StokesProblem class implementation} + // @sect3{StokesProblem class implementation} - // @sect4{StokesProblem::StokesProblem} - - // The constructor of this class looks very - // similar to the one of step-20. The constructor - // initializes the variables for the polynomial - // degree, triangulation, finite element system - // and the dof handler. The underlying polynomial - // functions are of order degree+1 for - // the vector-valued velocity components and - // of order degree in pressure. - // This gives the LBB-stable element pair - // Q(degree+1)Q(degree). - - // Note that we initialize the triangulation - // with a MeshSmoothing argument, which ensures - // that the refinement of cells is done - // in a way that the approximation of the - // PDE solution remains well-behaved (problems - // arise if grids are too unstructered), - // see the documentation of - // Triangulation::MeshSmoothing - // for details. + // @sect4{StokesProblem::StokesProblem} + + // The constructor of this class looks very + // similar to the one of step-20. The + // constructor initializes the variables for + // the polynomial degree, triangulation, + // finite element system and the dof + // handler. The underlying polynomial + // functions are of order + // degree+1 for the + // vector-valued velocity components and of + // order degree in pressure. + // This gives the LBB-stable element pair + // Q(degree+1)Q(degree). + // + // Note that we initialize the triangulation + // with a MeshSmoothing argument, which + // ensures that the refinement of cells is + // done in a way that the approximation of + // the PDE solution remains well-behaved + // (problems arise if grids are too + // unstructered), see the documentation of + // Triangulation::MeshSmoothing + // for details. template StokesProblem::StokesProblem (const unsigned int degree) : @@ -472,72 +468,72 @@ StokesProblem::StokesProblem (const unsigned int degree) {} - // @sect4{StokesProblem::setup_dofs} + // @sect4{StokesProblem::setup_dofs} - // Given a mesh, this function associates - // the degrees of freedom with it and - // creates the corresponding matrices and - // vectors. + // Given a mesh, this function associates + // the degrees of freedom with it and + // creates the corresponding matrices and + // vectors. template void StokesProblem::setup_dofs () { - // Release preconditioner from - // previous steps since it - // will definitely not be needed - // any more after this point. + // Release preconditioner from + // previous steps since it + // will definitely not be needed + // any more after this point. A_preconditioner.reset (); dof_handler.distribute_dofs (fe); - // In order to make the ILU preconditioner - // (in 3D) to work efficiently, the dofs - // are renumbered using the Cuthill-McKee - // algorithm. Though, the block structure - // of velocity and pressure shall be as in - // step-20. This is done in two steps. First, - // all dofs are renumbered by - // DoFRenumbering::Cuthill_McKee, - // and then we renumber once again by - // components. Since - // DoFRenumbering::component_wise - // does not touch the renumbering within - // the individual blocks, the basic - // renumbering from Cuthill-McKee remains. + // In order to make the ILU preconditioner + // (in 3D) to work efficiently, the dofs + // are renumbered using the Cuthill-McKee + // algorithm. Though, the block structure + // of velocity and pressure shall be as in + // step-20. This is done in two + // steps. First, all dofs are renumbered by + // DoFRenumbering::Cuthill_McKee, + // and then we renumber once again by + // components. Since + // DoFRenumbering::component_wise + // does not touch the renumbering within + // the individual blocks, the basic + // renumbering from Cuthill-McKee remains. DoFRenumbering::Cuthill_McKee (dof_handler); - // There is one more change: There - // is no reason in creating dim - // blocks for the velocity components, - // so they can all be grouped in only - // one block. The vector - // block_component does precisely - // this: velocity values correspond to block - // 0, and pressure values will sit in block - // 1. + // There is one more change: There is no + // reason in creating dim + // blocks for the velocity components, so + // they can all be grouped in only one + // block. The vector + // block_component does + // precisely this: velocity values + // correspond to block 0, and pressure + // values will sit in block 1. std::vector block_component (dim+1,0); block_component[dim] = 1; DoFRenumbering::component_wise (dof_handler, block_component); - // Since we use adaptively refined grids - // the constraint matrix for hanging - // node constraints is generated from - // the dof handler. + // Since we use adaptively refined grids + // the constraint matrix for hanging node + // constraints is generated from the dof + // handler. hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints); hanging_node_constraints.close (); - // In analogy to step-20, we count - // the dofs in the individual components. - // We could do this in the same way as - // there, but we want to operate on - // the block structure we used already for - // the renumbering: The function - // DoFTools::count_dofs_per_block - // does the same as - // DoFTools::count_dofs_per_component, - // but now grouped as velocity and - // pressure block via block_component. + // In analogy to step-20, we count the dofs + // in the individual components. We could + // do this in the same way as there, but we + // want to operate on the block structure + // we used already for the renumbering: The + // function + // DoFTools::count_dofs_per_block + // does the same as + // DoFTools::count_dofs_per_component, + // but now grouped as velocity and pressure + // block via block_component. std::vector dofs_per_block (2); DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); const unsigned int n_u = dofs_per_block[0], @@ -551,40 +547,37 @@ void StokesProblem::setup_dofs () << " (" << n_u << '+' << n_p << ')' << std::endl; - // Release the memory previously attached - // to the system matrix and untie it - // from the old sparsity pattern prior to - // generating the current data structure. + // Release the memory previously attached + // to the system matrix and untie it from + // the old sparsity pattern prior to + // generating the current data structure. system_matrix.clear (); - // The next task is to allocate a - // sparsity pattern for the system - // matrix we will create. We could do - // this in the same way as in step-20, - // though, there is a major reason - // not to do so. In 3D, the function - // DoFTools::max_couplings_between_dofs - // yields a very large number for the - // coupling between the individual dofs, - // so that the memory initially provided - // for the creation of the sparsity pattern - // of the matrix is far too much - so - // much actually that it won't even fit - // into the physical memory of most - // systems already for moderately-sized 3D - // problems. See also the discussion in - // step-18. - // Instead, we use a temporary object of - // the class - // BlockCompressedSparsityPattern, - // which is a block version of the - // compressed sparsity patterns from - // step-11 and step-18. All this is done - // inside a new scope, which means that - // the memory of csp will be - // released once the information has been - // copied to - // sparsity_pattern. + // The next task is to allocate a sparsity + // pattern for the system matrix we will + // create. We could do this in the same way + // as in step-20, though, there is a major + // reason not to do so. In 3D, the function + // DoFTools::max_couplings_between_dofs + // yields a very large number for the + // coupling between the individual dofs, so + // that the memory initially provided for + // the creation of the sparsity pattern of + // the matrix is far too much - so much + // actually that it won't even fit into the + // physical memory of most systems already + // for moderately-sized 3D problems. See + // also the discussion in step-18. + // Instead, we use a temporary object of + // the class + // BlockCompressedSparsityPattern, which is + // a block version of the compressed + // sparsity patterns from step-11 and + // step-18. All this is done inside a new + // scope, which means that the memory of + // csp will be released once + // the information has been copied to + // sparsity_pattern. { BlockCompressedSparsityPattern csp; @@ -604,10 +597,10 @@ void StokesProblem::setup_dofs () std::ofstream out ("sparsity_pattern.gpl"); sparsity_pattern.block(0,0).print_gnuplot(out); - // Finally, the system matrix, - // solution and right hand side are - // created from the block - // structure as in step-20. + // Finally, the system matrix, + // solution and right hand side are + // created from the block + // structure as in step-20. system_matrix.reinit (sparsity_pattern); solution.reinit (2); @@ -622,16 +615,16 @@ void StokesProblem::setup_dofs () } - // @sect4{StokesProblem::assemble_system} + // @sect4{StokesProblem::assemble_system} - // The assembly process follows the - // discussion in step-20 and in the - // introduction. We use the well-known - // abbreviations for the data structures - // that hold the local matrix, right - // hand side, and global - // numbers of the degrees of freedom - // for the present cell. + // The assembly process follows the + // discussion in step-20 and in the + // introduction. We use the well-known + // abbreviations for the data structures + // that hold the local matrix, right + // hand side, and global + // numbers of the degrees of freedom + // for the present cell. template void StokesProblem::assemble_system () { @@ -662,13 +655,13 @@ void StokesProblem::assemble_system () std::vector local_dof_indices (dofs_per_cell); - // As usual, we create objects that - // hold the functions for the right - // hand side and Neumann boundary - // function, and, additionally, - // an array that holds the respective - // function values at the quadrature - // points. + // As usual, we create objects that + // hold the functions for the right + // hand side and Neumann boundary + // function, and, additionally, + // an array that holds the respective + // function values at the quadrature + // points. const PressureBoundaryValues pressure_boundary_values; std::vector boundary_values (n_face_q_points); @@ -684,11 +677,11 @@ void StokesProblem::assemble_system () const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); - // This starts the loop over all - // cells. With the abbreviations - // extract_u etc. - // introduced above, it is - // evident what is going on. + // This starts the loop over all + // cells. With the abbreviations + // extract_u etc. + // introduced above, it is + // evident what is going on. typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -711,25 +704,29 @@ void StokesProblem::assemble_system () const double div_phi_i_u = fe_values[velocities].divergence (i, q); const double phi_i_p = fe_values[pressure].value (i, q); - for (unsigned int j=0; j phi_j_grads_u = fe_values[velocities].symmetric_gradient (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); - // Note the way we write the - // contributions - // phi_i_p * phi_j_p , - // yielding a pressure mass matrix, - // into the same data structure as - // the terms for the actual - // Stokes system - in accordance with - // the description in the introduction. - // They won't be mixed up, since - // phi_i_p * phi_j_p - // is only non-zero when all the - // other terms vanish (and the other - // way around). + // Note the way we write + // the contributions + // phi_i_p * phi_j_p + // , yielding a + // pressure mass matrix, + // into the same data + // structure as the terms + // for the actual Stokes + // system - in accordance + // with the description in + // the introduction. They + // won't be mixed up, since + // phi_i_p * + // phi_j_p is only + // non-zero when all the + // other terms vanish (and + // the other way around). // // Note also that operator* // is overloaded for @@ -738,27 +735,27 @@ void StokesProblem::assemble_system () // product between the two // tensors in the first // line: - 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) - * fe_values.JxW(q); - - } - const unsigned int component_i = - fe.system_to_component_index(i).first; - local_rhs(i) += fe_values.shape_value(i,q) * - rhs_values[q](component_i) * - fe_values.JxW(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) + * fe_values.JxW(q); + + } + const unsigned int component_i = + fe.system_to_component_index(i).first; + local_rhs(i) += fe_values.shape_value(i,q) * + rhs_values[q](component_i) * + fe_values.JxW(q); + } + } - // Here we add the contributions from - // Neumann (pressure) boundary conditions - // at faces on the domain boundary that - // have the boundary flag "0", i.e. those - // that are not subject to Dirichlet - // conditions. + // Here we add the contributions from + // Neumann (pressure) boundary + // conditions at faces on the domain + // boundary that have the boundary flag + // "0", i.e. those that are not subject + // to Dirichlet conditions. for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) @@ -784,16 +781,15 @@ void StokesProblem::assemble_system () } } - // The final step is, as usual, - // the transfer of the local - // contributions to the global - // system matrix. This works - // also in the case of block - // vectors and matrices, and - // also the terms constituting - // the pressure mass matrix are - // written into the correct position - // without any further interaction. + // The final step is, as usual, the + // transfer of the local contributions + // to the global system matrix. This + // works also in the case of block + // vectors and matrices, and also the + // terms constituting the pressure mass + // matrix are written into the correct + // position without any further + // interaction. cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i::assemble_system () system_rhs(local_dof_indices[i]) += local_rhs(i); } - // After the addition of the local - // contributions, we have to - // condense the hanging node - // constraints and interpolate - // Dirichlet boundary conditions. - // Note that Dirichlet boundary - // conditions are only condensed - // in boundary points that are - // labeled with "1", indicating - // that Dirichlet data is to be - // set. - // There is one more thing, though. - // The function describing the - // Dirichlet conditions was - // defined for all components, both - // velocity and pressure. However, - // the Dirichlet conditions are to - // be set for the velocity only. - // To this end, we use a - // component_mask that - // filters away the pressure - // component, so that the condensation - // is performed on - // velocity dofs. + // After the addition of the local + // contributions, we have to condense the + // hanging node constraints and interpolate + // Dirichlet boundary conditions. Note + // that Dirichlet boundary conditions are + // only condensed in boundary points that + // are labeled with "1", indicating that + // Dirichlet data is to be set. There is + // one more thing, though. The function + // describing the Dirichlet conditions was + // defined for all components, both + // velocity and pressure. However, the + // Dirichlet conditions are to be set for + // the velocity only. To this end, we use + // a component_mask that + // filters away the pressure component, so + // that the condensation is performed on + // velocity dofs. hanging_node_constraints.condense (system_matrix); hanging_node_constraints.condense (system_rhs); @@ -849,48 +839,43 @@ void StokesProblem::assemble_system () system_rhs); } - // Before we're going to solve - // this linear system, we generate - // a preconditioner for the - // velocity-velocity matrix, - // i.e., block(0,0) - // in the system matrix. As mentioned - // above, this depends on the - // spatial dimension. Since this - // handled automatically by the - // template - // in InnerPreconditioner, - // we don't have to manually - // intervene at this point any - // further. + // Before we're going to solve this linear + // system, we generate a preconditioner for + // the velocity-velocity matrix, i.e., + // block(0,0) in the system + // matrix. As mentioned above, this depends + // on the spatial dimension. Since this + // handled automatically by the template + // in + // InnerPreconditioner, we + // don't have to manually intervene at this + // point any further. std::cout << " Computing preconditioner..." << std::endl << std::flush; A_preconditioner = boost::shared_ptr::type>(new typename InnerPreconditioner::type()); A_preconditioner->initialize (system_matrix.block(0,0), - typename InnerPreconditioner::type::AdditionalData()); + typename InnerPreconditioner::type::AdditionalData()); } - // @sect4{StokesProblem::solve} + // @sect4{StokesProblem::solve} - // After the discussion in the - // introduction and the definition - // of the respective classes above, - // the implementation of the - // solve function is - // rather straigt-forward and done in - // a similar way as in step-20. To - // start with, we need an object of - // the InverseMatrix class - // that represents the inverse of - // the matrix A. As described in - // the introduction, the inverse - // is generated with the help - // of an inner preconditioner of - // type InnerPreconditioner. + // After the discussion in the introduction + // and the definition of the respective + // classes above, the implementation of the + // solve function is rather + // straigt-forward and done in a similar way + // as in step-20. To start with, we need an + // object of the InverseMatrix + // class that represents the inverse of the + // matrix A. As described in the + // introduction, the inverse is generated + // with the help of an inner preconditioner + // of type + // InnerPreconditioner. template void StokesProblem::solve () { @@ -899,16 +884,14 @@ void StokesProblem::solve () A_inverse (system_matrix.block(0,0), *A_preconditioner); Vector tmp (solution.block(0).size()); - // This is as in step-20. We generate - // the right hand side - // B A^{-1} F Ð G for the - // Schur complement and an object - // that represents the respective - // linear operation B A^{-1} B^T, - // now with a template parameter - // indicating the preconditioner - - // in accordance with the definition - // of the class. + // This is as in step-20. We generate the + // right hand side $B A^{-1} F - G$ for the + // Schur complement and an object that + // represents the respective linear + // operation $B A^{-1} B^T$, now with a + // template parameter indicating the + // preconditioner - in accordance with the + // definition of the class. { Vector schur_rhs (solution.block(1).size()); A_inverse.vmult (tmp, system_rhs.block(0)); @@ -918,90 +901,79 @@ void StokesProblem::solve () SchurComplement::type> schur_complement (system_matrix, A_inverse); - // The usual control structures for - // the solver call are created... + // The usual control structures for + // the solver call are created... SolverControl solver_control (system_matrix.block(0,0).m(), - 1e-6*schur_rhs.l2_norm()); + 1e-6*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - // Now to the preconditioner to the - // Schur complement. As explained in the - // introduction, the preconditioning - // is done by a mass matrix in the - // pressure variable. - // It is stored in the (1,1) block - // of the system matrix (that is not - // used elsewhere in this function). - - // Actually, the solver needs to have - // the preconditioner in the form - // P^{-1}, so we need to create - // an inverse operation. Once again, - // we use an object of the class - // InverseMatrix, which - // implements the vmult - // operation that is needed by the solver. - // In this case, we have to invert - // the pressure mass matrix. As it - // already turned out in earlier tutorial - // programs, the inversion of a mass - // matrix is a rather cheap and - // straight-forward operation (compared - // to, e.g., a Laplace matrix). The CG - // method with SSOR preconditioning - // converges in 10-20 steps, - // independently on the mesh size. - // This is precisely what we do here: - // We choose an SSOR preconditioner - // with parameter 1.2 and take it along - // to the InverseMatrix object via - // the corresponding template parameter. - // A CG solver is then called within - // the vmult operation of the inverse - // matrix. + // Now to the preconditioner to the Schur + // complement. As explained in the + // introduction, the preconditioning is + // done by a mass matrix in the pressure + // variable. It is stored in the $(1,1)$ + // block of the system matrix (that is + // not used elsewhere in this function). + // + // Actually, the solver needs to have the + // preconditioner in the form $P^{-1}$, so + // we need to create an inverse + // operation. Once again, we use an + // object of the class + // InverseMatrix, which + // implements the vmult + // operation that is needed by the + // solver. In this case, we have to + // invert the pressure mass matrix. As it + // already turned out in earlier tutorial + // programs, the inversion of a mass + // matrix is a rather cheap and + // straight-forward operation (compared + // to, e.g., a Laplace matrix). The CG + // method with SSOR preconditioning + // converges in 10-20 steps, + // independently on the mesh size. This + // is precisely what we do here: We + // choose an SSOR preconditioner with + // parameter 1.2 and take it along to the + // InverseMatrix object via the + // corresponding template parameter. A + // CG solver is then called within the + // vmult operation of the inverse matrix. PreconditionSSOR<> preconditioner; preconditioner.initialize (system_matrix.block(1,1), 1.2); InverseMatrix,PreconditionSSOR<> > m_inverse (system_matrix.block(1,1), preconditioner); - // With the Schur complement and an - // efficient preconditioner at hand, - // we can solve the respective - // equation in the usual way. - try - { - cg.solve (schur_complement, solution.block(1), schur_rhs, - m_inverse); - } - catch (...) - { - abort (); - } + // With the Schur complement and an + // efficient preconditioner at hand, + // we can solve the respective + // equation in the usual way. + cg.solve (schur_complement, solution.block(1), schur_rhs, + m_inverse); - // After this first solution step, - // the hanging node constraints have - // to be distributed to the solution - - // in order to achieve a consistent - // pressure field. + // After this first solution step, + // the hanging node constraints have + // to be distributed to the solution - + // in order to achieve a consistent + // pressure field. hanging_node_constraints.distribute (solution); std::cout << " " - << solver_control.last_step() - << " outer CG Schur complement iterations for pressure" - << std::flush - << std::endl; + << solver_control.last_step() + << " outer CG Schur complement iterations for pressure" + << std::flush + << std::endl; } - // As in step-20, we finally need to - // solve for the velocity equation - // where we plug in the the solution - // to the pressure equation. This involves - // only objects we already know - so - // we simply - // multiply p by B^T, subtract the - // right hand side and multiply - // by the inverse of A. + // As in step-20, we finally need to solve + // for the velocity equation where we plug + // in the the solution to the pressure + // equation. This involves only objects we + // already know - so we simply multiply p + // by $B^T$, subtract the right hand side and + // multiply by the inverse of A. { system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; @@ -1009,34 +981,30 @@ void StokesProblem::solve () A_inverse.vmult (solution.block(0), tmp); - // Again, we need to distribute - // the constraints from hanging nodes - // in order to obtain a constistent - // flow field. + // Again, we need to distribute the + // constraints from hanging nodes in + // order to obtain a constistent flow + // field. hanging_node_constraints.distribute (solution); } } - // @sect4{StokesProblem::output_results} + // @sect4{StokesProblem::output_results} - // The next function generates graphical - // output. In this example, we are going - // to use the VTK file format. - // We attach names to the individual - // variables in the problem - - // velocity to the dim - // components of velocity and - // p to the pressure. - // In order to tell the VTK file - // which components are vectors - // and which scalars, we need to - // add that information as well - - // achieved by the - // DataComponentInterpretation - // class. - // The rest of the function is - // then the same as in step-20. + // The next function generates graphical + // output. In this example, we are going to + // use the VTK file format. We attach names + // to the individual variables in the problem + // - velocity to the dim + // components of velocity and p + // to the pressure. In order to tell the VTK + // file which components are vectors and + // which scalars, we need to add that + // information as well - achieved by the + // DataComponentInterpretation + // class. The rest of the function is then + // the same as in step-20. template void StokesProblem::output_results (const unsigned int refinement_cycle) const @@ -1056,8 +1024,8 @@ StokesProblem::output_results (const unsigned int refinement_cycle) const = DataComponentInterpretation::component_is_part_of_vector; data_out.add_data_vector (solution, solution_names, - DataOut::type_dof_data, - data_component_interpretation); + DataOut::type_dof_data, + data_component_interpretation); data_out.build_patches (); @@ -1071,21 +1039,21 @@ StokesProblem::output_results (const unsigned int refinement_cycle) const } - // @sect4{StokesProblem::refine_mesh} + // @sect4{StokesProblem::refine_mesh} - // This is the last interesting function - // of the StokesProblem class. - // As indicated by its name, it takes the - // solution to the problem and - // refines the mesh where this is - // needed. The procedure is the same - // as in the respective step in - // step-6, with the exception that - // we base the refinement only on the - // change in pressure, i.e., we call - // the Kelly error estimator with a - // mask object. Additionally, we do - // not coarsen the grid again. + // This is the last interesting function + // of the StokesProblem class. + // As indicated by its name, it takes the + // solution to the problem and + // refines the mesh where this is + // needed. The procedure is the same + // as in the respective step in + // step-6, with the exception that + // we base the refinement only on the + // change in pressure, i.e., we call + // the Kelly error estimator with a + // mask object. Additionally, we do + // not coarsen the grid again. template void StokesProblem::refine_mesh () @@ -1108,23 +1076,25 @@ StokesProblem::refine_mesh () } - // @sect4{StokesProblem::run} + // @sect4{StokesProblem::run} - // The last step in the Stokes class - // is, as usual, the program that generates - // the initial grid and calls the other - // functions in the respective order. + // The last step in the Stokes class + // is, as usual, the program that generates + // the initial grid and calls the other + // functions in the respective order. template void StokesProblem::run () { - // We start off with a rectangle of - // size 4 x 1 (x 1), placed in R^2/R^3 - // as (-2,2)x(-1,0) or (-2,2)x(0,1)x(-1,1), - // respectively. It is natural to start - // with equal mesh size in each direction, - // so we subdivide the initial rectangle - // four times in the first coordinate - // direction. + // We start off with a rectangle of size $4 + // \times 1$ (in 2d) or $4 \times 1 times + // 1$ (in 3d), placed in $R^2/R^3$ as + // $(-2,2)times(-1,0)$ or + // $(-2,2)\times(0,1)\times(-1,1)$, + // respectively. It is natural to start + // with equal mesh size in each direction, + // so we subdivide the initial rectangle + // four times in the first coordinate + // direction. std::vector subdivisions (dim, 1); subdivisions[0] = 4; @@ -1137,38 +1107,37 @@ void StokesProblem::run () Point(2,0) : Point(2,1,0))); - // A boundary indicator is set to all - // boundaries that are subject to - // Dirichlet boundary conditions, i.e. - // to faces that are located at 0 in - // the last coordinate direction. See - // the example description above for - // details. + // A boundary indicator is set to all + // boundaries that are subject to Dirichlet + // boundary conditions, i.e. to faces that + // are located at 0 in the last coordinate + // direction. See the example description + // above for details. for (typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(); + cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) for (unsigned int f=0; f::faces_per_cell; ++f) if (cell->face(f)->center()[dim-1] == 0) - { - cell->face(f)->set_boundary_indicator(1); - /* - for (unsigned int e=0; e::lines_per_face; ++e) - cell->face(f)->line(e)->set_boundary_indicator (1); - */ - } + { + cell->face(f)->set_boundary_indicator(1); + /* + for (unsigned int e=0; e::lines_per_face; ++e) + cell->face(f)->line(e)->set_boundary_indicator (1); + */ + } - // We employ an initial refinement before - // solving for the first time. In 3D, - // there are going to be more dofs, so - // we refine less there. + // We employ an initial refinement before + // solving for the first time. In 3D, there + // are going to be more dofs, so we refine + // less there. triangulation.refine_global (4-dim); - // As first seen in step-6, we cycle - // over the different refinement levels - // and refine (if not the first step), - // setup the dofs and matrices, assemble, - // solve and create an output. + // As first seen in step-6, we cycle over + // the different refinement levels and + // refine (if not the first step), setup + // the dofs and matrices, assemble, solve + // and create an output. for (unsigned int refinement_cycle = 0; refinement_cycle<7; ++refinement_cycle) { @@ -1192,13 +1161,12 @@ void StokesProblem::run () } - // @sect3{The main function} + // @sect3{The main function} - // The main function is the same as - // in step-20. We pass the element - // degree as a parameter and - // choose the space dimension at the - // well-known template slot. + // The main function is the same as in + // step-20. We pass the element degree as a + // parameter and choose the space dimension + // at the well-known template slot. int main () { try -- 2.39.5