From 8aa9e95fef5ac076fac086aa642f67c15ca29948 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 26 Feb 2008 16:04:50 +0000 Subject: [PATCH] Go over the comment in the text, work on them and simplify the code in a few places. Also try to keep more code and comments together to produce slightly bigger code blocks that make grasping the flow of control a bit simpler. git-svn-id: https://svn.dealii.org/trunk@15790 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 711 ++++++++++++++++------------ 1 file changed, 409 insertions(+), 302 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 855329f5a9..0310047b70 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -165,28 +165,40 @@ class StokesProblem // @sect3{%Boundary values and right hand side} - // As in step-20 and most other example - // programs, the next task is to define the - // data for the PDE: For the - // Stokes problem, we are going to use - // natural 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. + // As in step-20 and most other + // example programs, the next task is + // to define the data for the PDE: + // For the Stokes problem, we are + // going to use natural boundary + // values on parts of the boundary + // (i.e. homogenous Neumann-type) for + // which we won't have to do anything + // special (the homogeneity implies + // that the corresponding terms in + // the weak form are simply zero), + // and boundary conditions on the + // velocity (Dirichlet-type) on the + // rest of the boundary, as described + // in the introduction. + // + // In order to enforce the Dirichlet + // boundary values on the velocity, + // we will use the + // VectorTools::interpolate_boundary_values + // function as usual which requires + // us to write a function object with + // as many components as the finite + // element has. In other words, we + // have to define the function on the + // $(u,p)$-space, but we are going to + // filter out the pressure component + // when interpolating the boundary + // values. - // Given the problem described in the - // introduction, we know which values to - // set for the respective functions. + // The following function object is a + // representation of the boundary + // values described in the + // introduction: template class BoundaryValues : public Function { @@ -206,6 +218,9 @@ double BoundaryValues::value (const Point &p, const unsigned int component) const { + Assert (component < this->n_components, + ExcIndexRange (component, 0, this->n_components)); + if (component == 0) return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0)); return 0; @@ -223,8 +238,9 @@ BoundaryValues::vector_value (const Point &p, - // We implement similar functions - // for the right hand side. + // We implement similar functions for + // the right hand side which for the + // current example is simply zero: template class RightHandSide : public Function { @@ -268,20 +284,31 @@ RightHandSide::vector_value (const Point &p, // @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. + // The InverseMatrix + // class represents the data + // structure for an inverse + // matrix. It is derived from the one + // in step-20. The only difference is + // that we now do include a + // preconditioner to the matrix since + // we will apply this class to + // different kinds of matrices that + // will require different + // preconditioners (in step-20 we did + // not use a preconditioner in this + // class at all). The types of matrix + // and preconditioner are passed to + // this class via template + // parameters, and matrix and + // preconditioner objects of these + // types will then be passed to the + // constructor 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 { @@ -294,9 +321,7 @@ class InverseMatrix : public Subscriptor private: const SmartPointer matrix; - const Preconditioner &preconditioner; - - mutable GrowingVectorMemory<> vector_memory; + const SmartPointer preconditioner; }; @@ -305,7 +330,7 @@ InverseMatrix::InverseMatrix (const Matrix &m, const Preconditioner &preconditioner) : matrix (&m), - preconditioner (preconditioner) + preconditioner (&preconditioner) {} @@ -313,53 +338,32 @@ InverseMatrix::InverseMatrix (const Matrix &m, // 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. - - // 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. + // Note that 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. template void InverseMatrix::vmult (Vector &dst, const Vector &src) const { SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); - SolverCG<> cg (solver_control, vector_memory); + SolverCG<> cg (solver_control); dst = 0; - cg.solve (*matrix, dst, src, preconditioner); + cg.solve (*matrix, dst, src, *preconditioner); } @@ -375,22 +379,22 @@ void InverseMatrix::vmult (Vector &dst, // consequence of the definition above, the // declaration InverseMatrix now // contains the second template parameter - // from preconditioning as above, which - // effects the SmartPointer + // for a preconditioner class as above, which + // affects the SmartPointer // object m_inverse as well. template class SchurComplement : public Subscriptor { public: - SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix, Preconditioner> &Minv); + SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix, Preconditioner> &A_inverse); void vmult (Vector &dst, const Vector &src) const; private: const SmartPointer > system_matrix; - const SmartPointer, Preconditioner> > m_inverse; + const SmartPointer, Preconditioner> > A_inverse; mutable Vector tmp1, tmp2; }; @@ -399,13 +403,13 @@ class SchurComplement : public Subscriptor template SchurComplement:: -SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix,Preconditioner> &Minv) +SchurComplement (const BlockSparseMatrix &system_matrix, + const InverseMatrix,Preconditioner> &A_inverse) : - system_matrix (&A), - m_inverse (&Minv), - tmp1 (A.block(0,0).m()), - tmp2 (A.block(0,0).m()) + system_matrix (&system_matrix), + A_inverse (&A_inverse), + tmp1 (system_matrix.block(0,0).m()), + tmp2 (system_matrix.block(0,0).m()) {} @@ -414,7 +418,7 @@ void SchurComplement::vmult (Vector &dst, const Vector &src) const { system_matrix->block(0,1).vmult (tmp1, src); - m_inverse->vmult (tmp2, tmp1); + A_inverse->vmult (tmp2, tmp1); system_matrix->block(1,0).vmult (dst, tmp2); } @@ -423,18 +427,22 @@ void SchurComplement::vmult (Vector &dst, // @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, + // 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). + // vector-valued velocity components + // and of order degree + // for the pressure. This gives the + // LBB-stable element pair + // $Q_{degree+1}^d\times Q_{degree}$, + // often referred to as the + // Taylor-Hood element. // // Note that we initialize the triangulation // with a MeshSmoothing argument, which @@ -458,53 +466,88 @@ StokesProblem::StokesProblem (const unsigned int degree) // @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. At the beginning it also + // releases the pointer to the + // preconditioner object (if the + // shared pointer pointed at anything + // at all at this point) since it + // will definitely not be needed any + // more after this point and will + // have to be re-computed after + // assembling the matrix, and unties + // the sparse matrix from its + // sparsity pattern object. + // + // We the procedd with distributing + // degrees of freedom and renumbering + // them: In order to make the ILU + // preconditioner (in 3D) work + // efficiently, the degrees of + // freedom are renumbered using the + // Cuthill-McKee algorithm as this + // reduces the bandwidth of the + // matrix. On the other hand, we need + // to preserve the block structure of + // velocity and pressure already seen + // in in step-20 and step-21. 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. + // + // There is one more change compared + // to previous tutorial programs: + // There is no reason in sorting the + // dim velocity + // components individually. In fact, + // rather than first enumerating all + // $x$-velocities, then all + // $y$-velocities, etc, we would like + // to keep all velocities at the same + // location together and only + // separate between velocities (all + // components) and pressures. By + // default, this is not what the + // DoFRenumbering::component_wise + // function does: it treats each + // vector component separately; what + // we have to do is group several + // components into "blocks" and pass + // this block structure to that + // function. Consequently, we + // allocate a vector + // block_component with + // as many elements as there are + // components and describe all + // velocity components to correspond + // to block 0, while the pressure + // component will form block 1: template void StokesProblem::setup_dofs () { - // Release preconditioner from - // previous steps since it - // will definitely not be needed - // any more after this point. A_preconditioner.reset (); + system_matrix.clear (); - 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. + dof_handler.distribute_dofs (fe); 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. 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 + // constraints is generated from the DoF // handler. hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, @@ -534,42 +577,43 @@ void StokesProblem::setup_dofs () << dof_handler.n_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. - 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 + // 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 + // yields a conservative, 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 + // the initial sparsity pattern + // 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; + BlockCompressedSparsityPattern csp (2,2); - csp.reinit (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(1,0).reinit (n_p, n_u); csp.block(0,1).reinit (n_u, n_p); @@ -585,7 +629,7 @@ void StokesProblem::setup_dofs () // Finally, the system matrix, // solution and right hand side are // created from the block - // structure as in step-20. + // structure as in step-20: system_matrix.reinit (sparsity_pattern); solution.reinit (2); @@ -650,11 +694,8 @@ 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. + // We can then start the loop over all + // cells: typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -683,31 +724,7 @@ void StokesProblem::assemble_system () 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 also that operator* - // is overloaded for - // symmetric tensors, - // yielding the scalar - // 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 @@ -721,7 +738,33 @@ void StokesProblem::assemble_system () rhs_values[q](component_i) * fe_values.JxW(q); } - } + } + + // Note that in the above + // computation of the local + // matrix contribution we added + // the term phi_i_p * + // phi_j_p , yielding a + // pressure mass matrix in the + // $(1,1)$ block of the matrix + // as discussed in the + // introduction. That this term + // only ends up in the $(1,1)$ + // block stems from the fact + // that both of the factors in + // phi_i_p * + // phi_j_p are only + // non-zero when all the other + // terms vanish (and the other + // way around). + // + // Note also that operator* is + // overloaded for symmetric + // tensors, yielding the scalar + // product between the two + // tensors in the first line of + // the local matrix + // contribution. // The final step is, as usual, the // transfer of the local contributions @@ -731,7 +774,7 @@ void StokesProblem::assemble_system () // terms constituting the pressure mass // matrix are written into the correct // position without any further - // interaction. + // interaction: cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i::assemble_system () } // 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. + // contributions, we have to + // condense the hanging node + // constraints and interpolate + // Dirichlet boundary conditions. + // Further down below where we set + // up the mesh, we will associate + // the top boundary where we impose + // Dirichlet boundary conditions + // with boundary indicator 1. We + // will have to pass this boundary + // indicator as second argument to + // the function below interpolating + // boundary values. 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 out the pressure + // component, so that the + // condensation is performed on + // velocity degrees of freedom + // only: hanging_node_constraints.condense (system_matrix); hanging_node_constraints.condense (system_rhs); @@ -781,17 +834,21 @@ 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 - // dim 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 the two + // classes described by the + // InnerPreconditioner@ :: type + // typedef have the same interface, + // we do not have to do anything + // different whether we want to use + // a sparse direct solver or an + // ILU: std::cout << " Computing preconditioner..." << std::endl << std::flush; A_preconditioner @@ -817,7 +874,7 @@ void StokesProblem::assemble_system () // introduction, the inverse is generated // with the help of an inner preconditioner // of type - // InnerPreconditioner. + // InnerPreconditioner@::type. template void StokesProblem::solve () { @@ -849,13 +906,16 @@ void StokesProblem::solve () 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). + // 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 + // anywhere else but in + // preconditioning). // // Actually, the solver needs to have the // preconditioner in the form $P^{-1}$, so @@ -888,16 +948,19 @@ void StokesProblem::solve () 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. + // With the Schur complement and + // an efficient preconditioner at + // hand, we can solve the + // respective equation for the + // pressure (i.e. block 0 in the + // solution vector) 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 - + // to be distributed to the solution // in order to achieve a consistent // pressure field. hanging_node_constraints.distribute (solution); @@ -909,24 +972,26 @@ void StokesProblem::solve () << 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 + // 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$. At the end, we need to + // distribute the constraints from + // hanging nodes in order to obtain + // a constistent flow field: { system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; tmp += system_rhs.block(0); 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. + hanging_node_constraints.distribute (solution); } } @@ -937,38 +1002,63 @@ void StokesProblem::solve () // 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 + // to the individual variables in the problem: + // velocity to the dim + // components of velocity and pressure + // to the pressure. + // + // Not all visualization programs + // have the ability to group + // individual vector components into + // a vector to provide vector plots; + // in particular, this holds for some + // VTK-based visualization + // programs. In this case, the + // logical grouping of components + // into vectors should already be + // described in the file containing + // the data. In other words, what we + // need to do is provide our output + // writers with a way to know which + // of the components of the finite + // element logically form a vector + // (with $d$ components in $d$ space + // dimensions) rather than letting + // them assume that we simply have a + // bunch of scalar fields. This is + // achieved using the members of the // DataComponentInterpretation - // class. The rest of the function is then + // namespace: as with the filename, + // we create a vector in which the + // first dim components + // refer to the velocities and are + // given the tag + // DataComponentInterpretation::component_is_part_of_vector; + // we finally push one tag + // DataComponentInterpretation::component_is_scalar + // to describe the grouping of the + // pressure variable. + + // The rest of the function is then // the same as in step-20. template void StokesProblem::output_results (const unsigned int refinement_cycle) const { std::vector solution_names (dim, "velocity"); - solution_names.push_back ("p"); + solution_names.push_back ("pressure"); - DataOut data_out; - - data_out.attach_dof_handler (dof_handler); - std::vector data_component_interpretation - (dim+1, DataComponentInterpretation::component_is_scalar); - for (unsigned int i=0; i data_out; + data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, solution_names, DataOut::type_dof_data, data_component_interpretation); - data_out.build_patches (); std::ostringstream filename; @@ -995,7 +1085,7 @@ StokesProblem::output_results (const unsigned int refinement_cycle) const // change in pressure, i.e., we call // the Kelly error estimator with a // mask object. Additionally, we do - // not coarsen the grid again. + // not coarsen the grid again: template void StokesProblem::refine_mesh () @@ -1021,35 +1111,48 @@ StokesProblem::refine_mesh () // @sect4{StokesProblem::run} // The last step in the Stokes class - // is, as usual, the program that generates + // is, as usual, the function that generates // the initial grid and calls the other // functions in the respective order. + // + // 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. To + // limit the scope of the variables + // involved in the creation of the + // mesh to the range where we + // actually need them, we put the + // entire block between a pair of + // braces: template void StokesProblem::run () { - // 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; + { + std::vector subdivisions (dim, 1); + subdivisions[0] = 4; + + const Point bottom_left = (dim == 2 ? + Point(-2,-1) : + Point(-2,0,-1)); + const Point top_right = (dim == 2 ? + Point(2,0) : + Point(2,1,0)); - GridGenerator::subdivided_hyper_rectangle (triangulation, - subdivisions, - (dim == 2 ? - Point(-2,-1) : - Point(-2,0,-1)), - (dim == 2 ? - Point(2,0) : - Point(2,1,0))); + GridGenerator::subdivided_hyper_rectangle (triangulation, + subdivisions, + bottom_left, + top_right); + } - // A boundary indicator is set to all + // A boundary indicator of 1 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 @@ -1069,17 +1172,21 @@ void StokesProblem::run () } - // 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 then apply an initial + // refinement before solving for + // the first time. In 3D, there are + // going to be more degrees of + // freedom, 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 + // (except for the first cycle), + // setup the degrees of freedom and + // matrices, assemble, solve and + // create output: for (unsigned int refinement_cycle = 0; refinement_cycle<7; ++refinement_cycle) { -- 2.39.5