From 33c2e3d158a40484454e4242f485f162b75e4acb Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 20 Feb 2008 14:37:26 +0000 Subject: [PATCH] Added comments to step-31.cc git-svn-id: https://svn.dealii.org/trunk@15746 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 857 ++++++++++++++++++++++------ 1 file changed, 689 insertions(+), 168 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index c4f08839d3..e7ea853b3d 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -1,7 +1,7 @@ -/* $Id: step-22.cc 15679 2008-01-24 23:28:37Z bangerth $ */ +/* $Id: step-31.cc 15679 2008-01-24 23:28:37Z bangerth $ */ /* Author: Wolfgang Bangerth, Texas A&M University, 2008 */ -/* $Id: step-22.cc 15679 2008-01-24 23:28:37Z bangerth $ */ +/* $Id: step-31.cc 15679 2008-01-24 23:28:37Z bangerth $ */ /* Version: $Name$ */ /* */ /* Copyright (C) 2008 by the deal.II authors */ @@ -12,7 +12,10 @@ /* further information on this license. */ - + // @sect3{Include files} + + // As usual, we start by including + // some well-known files. #include #include #include @@ -23,8 +26,6 @@ #include #include #include -#include -#include #include #include @@ -41,34 +42,62 @@ #include #include -#include #include #include #include -#include #include #include #include #include -#include + // 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. +#include + + // This is C++: #include #include + // As in all programs, the namespace is set + // to dealii. using namespace dealii; - + // @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. template struct InnerPreconditioner; + // 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@. template <> struct InnerPreconditioner<3> { @@ -76,7 +105,14 @@ struct InnerPreconditioner<3> }; - + // @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. template class StokesProblem { @@ -105,13 +141,38 @@ 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. 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. + + // Given the problem described in the + // introduction, we know which values to + // set for the respective functions. template class PressureBoundaryValues : public Function { @@ -169,6 +230,8 @@ BoundaryValues::vector_value (const Point &p, + // We implement similar functions + // for the right hand side. template class RightHandSide : public Function { @@ -180,6 +243,7 @@ class RightHandSide : public Function virtual void vector_value (const Point &p, Vector &value) const; + }; @@ -203,9 +267,22 @@ RightHandSide::vector_value (const Point &p, + // @sect3{extract_u and friends} - - + // The next four functions are needed for + // the assembly of the system matrix and + // the right hand side. They are very similar + // to the ones used in step-20, except + // that we are going to use Q(p+1)Qp elements + // instead of divergence-free Raviart-Thomas + // elements, which simplifies this procedure. + // The only function that is new is + // extract_grad_s_u, which + // gets the symmetric gradient of u. + // As discussed in the introduction, this + // is a second-rank tensor, formed by + // contributions from the gradient and its + // transpose. template Tensor<1,dim> extract_u (const FEValuesBase &fe_values, @@ -241,7 +318,6 @@ extract_div_u (const FEValuesBase &fe_values, } - template Tensor<2,dim> extract_grad_s_u (const FEValuesBase &fe_values, @@ -282,17 +358,58 @@ double extract_p (const FEValuesBase &fe_values, return 0; } + // @sect4{Inner product of second-rank tensors} + + // In the assembly process, we will need + // to form inner products of second-rank + // tensors. The way how to do this was + // discussed in the introduction - just + // take the sum of the product of the + // individual entries. +template +double +scalar_product (const Tensor<2,dim> &a, + const Tensor<2,dim> &b) +{ + double tmp = 0; + for (unsigned int i=0; iInverseMatrix 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. template class InverseMatrix : public Subscriptor { public: InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner); + const Preconditioner &preconditioner); void vmult (Vector &dst, - const Vector &src) const; + const Vector &src) const; private: const SmartPointer matrix; @@ -304,17 +421,62 @@ 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: + + // 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. 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); @@ -332,16 +494,31 @@ 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. 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; @@ -355,18 +532,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); @@ -374,52 +551,152 @@ void SchurComplement::vmult (Vector &dst, } - + // @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 discussion of + // Triangulation::MeshSmoothing + // for details. template StokesProblem::StokesProblem (const unsigned int degree) : degree (degree), - triangulation (Triangulation::maximum_smoothing), + triangulation (Triangulation::maximum_smoothing), fe (FE_Q(degree+1), dim, FE_Q(degree), 1), dof_handler (triangulation) {} - - + // @sect4{StokesProblem::setup_dofs} + + // 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 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. DoFRenumbering::Cuthill_McKee (dof_handler); - DoFRenumbering::component_wise (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 handler. hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, - hanging_node_constraints); + hanging_node_constraints); hanging_node_constraints.close (); - std::vector dofs_per_component (dim+1); - 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]; + // 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], + n_p = dofs_per_block[1]; std::cout << " Number of active cells: " << triangulation.n_active_cells() << std::endl << " Number of degrees of freedom: " << dof_handler.n_dofs() - << " (" << n_u << '+' << n_p <<')' + << " (" << n_u << '+' << n_p << ')' << std::endl; - + + // Clear the system matrix prior to + // generating the entries. 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 + // max_couplings_between_dofs@ + // yields a very large number for the + // coupling between the individual dofs, + // so that the memory initially provided + // in the reinit 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 discussing 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; @@ -436,6 +713,10 @@ void StokesProblem::setup_dofs () sparsity_pattern.copy_from (csp); } + // 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); @@ -450,19 +731,16 @@ void StokesProblem::setup_dofs () } -template -double -scalar_product (const Tensor<2,dim> &a, - const Tensor<2,dim> &b) -{ - double tmp = 0; - for (unsigned int i=0; i void StokesProblem::assemble_system () @@ -474,13 +752,15 @@ void StokesProblem::assemble_system () QGauss face_quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - update_gradients); + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); FEFaceValues fe_face_values (fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | update_JxW_values); + update_values | + update_normal_vectors | + update_quadrature_points | + update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; @@ -492,10 +772,26 @@ 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. const PressureBoundaryValues pressure_boundary_values; std::vector boundary_values (n_face_q_points); - + + const RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); + + // This starts the loop over all + // cells. With the abbreviations + // extract_u etc. + // introduced above, it is very + // clear what is going on. typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -504,38 +800,63 @@ void StokesProblem::assemble_system () fe_values.reinit (cell); local_matrix = 0; local_rhs = 0; - + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + for (unsigned int q=0; q phi_i_u = extract_u (fe_values, i, q); - + { + for (unsigned int i=0; i phi_i_grads_u= extract_grad_s_u (fe_values, i, q); const double div_phi_i_u = extract_div_u (fe_values, i, q); const double phi_i_p = extract_p (fe_values, i, q); + + for (unsigned int j=0; j phi_j_grads_u= extract_grad_s_u (fe_values, j, q); + const double div_phi_j_u = extract_div_u (fe_values, j, q); + const double phi_j_p = extract_p (fe_values, j, q); + - for (unsigned int j=0; j phi_j_grads_u = extract_grad_s_u (fe_values, j, q); - const double div_phi_j_u = extract_div_u (fe_values, j, q); - const double phi_j_p = extract_p (fe_values, j, q); - - local_matrix(i,j) += (scalar_product(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); - } + // Note how 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. + local_matrix(i,j) += (scalar_product(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. for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - if (cell->at_boundary(face_no)) + if (cell->at_boundary(face_no) && + (cell->face(face_no)->boundary_indicator() == 0)) { fe_face_values.reinit (cell, face_no); @@ -543,31 +864,65 @@ void StokesProblem::assemble_system () .value_list (fe_face_values.get_quadrature_points(), boundary_values); - for (unsigned int q=0; q - phi_i_u = extract_u (fe_face_values, i, q); - - local_rhs(i) += -(phi_i_u * - fe_face_values.normal_vector(q) * - boundary_values[q] * - fe_face_values.JxW(q)); - } + { + const Tensor<1,dim> + phi_i_u = extract_u (fe_face_values, i, q); + + local_rhs(i) += -(phi_i_u * + fe_face_values.normal_vector(q) * + boundary_values[q] * + fe_face_values.JxW(q)); + } } + // 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; icomponent_mask that + // filters away the pressure + // componenent, so that the condensation + // is performed only on + // velocity dofs. hanging_node_constraints.condense (system_matrix); hanging_node_constraints.condense (system_rhs); @@ -576,88 +931,203 @@ void StokesProblem::assemble_system () std::vector component_mask (dim+1, true); component_mask[dim] = false; VectorTools::interpolate_boundary_values (dof_handler, - 1, - BoundaryValues(), - boundary_values, - component_mask); + 1, + BoundaryValues(), + boundary_values, + component_mask); MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); + system_matrix, + solution, + system_rhs); } - std::cout << " Computing preconditioner..." << std::flush; + // 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()); - std::cout << std::endl; } + // @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. template void StokesProblem::solve () { - const InverseMatrix,typename InnerPreconditioner::type> + const InverseMatrix, + typename InnerPreconditioner::type> A_inverse (system_matrix.block(0,0), *A_preconditioner); Vector tmp (solution.block(0).size()); - Vector schur_rhs (solution.block(1).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. { + Vector schur_rhs (solution.block(1).size()); A_inverse.vmult (tmp, system_rhs.block(0)); system_matrix.block(1,0).vmult (schur_rhs, tmp); schur_rhs -= system_rhs.block(1); - + SchurComplement::type> schur_complement (system_matrix, A_inverse); + // 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()); SolverCG<> cg (solver_control); + // Now to the preconditioner to the + // Schur complement. As derived 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 + // program, 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 simple preconditioning + // with SSOR 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. 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); + cg.solve (schur_complement, solution.block(1), schur_rhs, + m_inverse); } catch (...) { - abort (); + abort (); } - - // produce a consistent flow field + + // After this first solution step, + // the hanging node constraints have + // to be distributed to the solution - + // that a consistent pressure field + // is achieved. hanging_node_constraints.distribute (solution); - std::cout << " " + std::cout << " " << solver_control.last_step() - << " CG Schur complement iterations for pressure." - << std::endl; + << " outer CG iterations for p " + << std::flush; } - + + // As in step-20, we finally need to + // solve for the velocity equation + // with the solution of the pressure + // equation at hand. We do not perform + // any direct solution of a linear + // system, but only need to + // 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; tmp += system_rhs.block(0); - + A_inverse.vmult (solution.block(0), tmp); - - // produce a consistent pressure 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} + + // 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 @@ -665,7 +1135,7 @@ StokesProblem::output_results (const unsigned int refinement_cycle) const std::vector solution_names (dim, "velocity"); solution_names.push_back ("p"); - DataOut data_out; + DataOut data_out; data_out.attach_dof_handler (dof_handler); @@ -677,22 +1147,36 @@ 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 (); std::ostringstream filename; filename << "solution-" - << Utilities::int_to_string (refinement_cycle, 2) - << ".vtk"; + << Utilities::int_to_string (refinement_cycle, 2) + << ".vtk"; std::ofstream output (filename.str().c_str()); data_out.write_vtk (output); } - + // @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. template void StokesProblem::refine_mesh () @@ -702,72 +1186,109 @@ StokesProblem::refine_mesh () std::vector component_mask (dim+1, false); component_mask[dim] = true; KellyErrorEstimator::estimate (dof_handler, - QGauss(degree+1), - typename FunctionMap::type(), - solution, - estimated_error_per_cell, - component_mask); + QGauss(degree+1), + typename FunctionMap::type(), + solution, + estimated_error_per_cell, + component_mask); GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.0); + estimated_error_per_cell, + 0.3, 0.0); triangulation.execute_coarsening_and_refinement (); } + // @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. 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. std::vector subdivisions (dim, 1); subdivisions[0] = 4; - + GridGenerator::subdivided_hyper_rectangle (triangulation, - subdivisions, - (dim == 2 ? - Point(-2,-1) : - Point(-2,0,-1)), - (dim == 2 ? - Point(2,0) : - Point(2,1,0))); + subdivisions, + (dim == 2 ? + Point(-2,-1) : + Point(-2,0,-1)), + (dim == 2 ? + 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. 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. 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. for (unsigned int refinement_cycle = 0; refinement_cycle<7; ++refinement_cycle) { std::cout << "Refinement cycle " << refinement_cycle << std::endl; if (refinement_cycle > 0) - refine_mesh (); + refine_mesh (); setup_dofs (); - std::cout << " Assembling..." << std::endl; + std::cout << " Assembling..." << std::endl << std::flush; assemble_system (); - std::cout << " Solving..." << std::endl; + std::cout << " Solving..." << std::flush; solve (); output_results (refinement_cycle); - std::cout << std::endl; + std::cout << std::endl << std::endl; } } - + // @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. int main () { try -- 2.39.5