/* 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 <base/quadrature_lib.h>
#include <base/logstream.h>
#include <base/function.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
- // 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 <lac/sparse_direct.h>
- // 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 <lac/sparse_ilu.h>
- // This is C++:
+ // This is C++:
#include <fstream>
#include <sstream>
- // 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 <int dim>
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 <code>SparseILU@<double></code>.
+ // And the ILU preconditioning in 3D, called
+ // by <code>SparseILU@<double></code>.
template <>
struct InnerPreconditioner<3>
{
};
- // @sect3{The <code>StokesProblem</code> class template}
+ // @sect3{The <code>StokesProblem</code> 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 <int dim>
class StokesProblem
{
BlockVector<double> solution;
BlockVector<double> 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<typename InnerPreconditioner<dim>::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 <code>assemble_system</code>.
+ // @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
+ // <code>assemble_system</code>.
- // 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 <int dim>
class PressureBoundaryValues : public Function<dim>
{
template <int dim>
double
BoundaryValues<dim>::value (const Point<dim> &p,
- const unsigned int component) const
+ const unsigned int component) const
{
if (component == 0)
return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
template <int dim>
void
BoundaryValues<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
+ Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = BoundaryValues<dim>::value (p, c);
- // We implement similar functions
- // for the right hand side.
+ // We implement similar functions
+ // for the right hand side.
template <int dim>
class RightHandSide : public Function<dim>
{
}
- // @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 <code>InverseMatrix</code> class template}
+ // @sect4{The <code>InverseMatrix</code> 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
- // <code>class Preconditioner</code>, so
- // the preconditioner type will be set when
- // an <code>InverseMatrix</code> object is
- // created. The member function
- // <code>vmult</code> 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 <code>class
+ // Preconditioner</code>, so the
+ // preconditioner type will be set when an
+ // <code>InverseMatrix</code> object is
+ // created. The member function
+ // <code>vmult</code> is, as in step-20, a
+ // multiplication with a vector, obtained by
+ // solving a linear system.
template <class Matrix, class Preconditioner>
class InverseMatrix : public Subscriptor
{
template <class Matrix, class Preconditioner>
InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
- const Preconditioner &preconditioner)
- :
- matrix (&m),
- preconditioner (preconditioner)
+ const Preconditioner &preconditioner)
+ :
+ matrix (&m),
+ preconditioner (preconditioner)
{}
- // This is the implementation of the
- // <code>vmult</code> function. We note
- // two things:
+ // This is the implementation of the
+ // <code>vmult</code> 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
- // <code>InverseMatrix</code> 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
- // <code>catch (SolverControl::NoConvergence) {}</code>
- // 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
+ // <code>InverseMatrix</code> 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 <code>catch
+ // (SolverControl::NoConvergence) {}</code>
+ // 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 <class Matrix, class Preconditioner>
void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double> &dst,
- const Vector<double> &src) const
+ const Vector<double> &src) const
{
SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
SolverCG<> cg (solver_control, vector_memory);
}
- // @sect4{The <code>SchurComplement</code> 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 <code>Preconditioner</code>
- // in order to access that when specifying
- // the respective type of the inverse
- // matrix class. As a consequence of the
- // definition above, the declaration
- // <code>InverseMatrix</code>
- // now contains the second template
- // parameter from preconditioning as above,
- // which effects the <code>SmartPointer</code>
- // object <code>m_inverse</code> as well.
+ // @sect4{The <code>SchurComplement</code> 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
+ // <code>Preconditioner</code> in order to
+ // access that when specifying the respective
+ // type of the inverse matrix class. As a
+ // consequence of the definition above, the
+ // declaration <code>InverseMatrix</code> now
+ // contains the second template parameter
+ // from preconditioning as above, which
+ // effects the <code>SmartPointer</code>
+ // object <code>m_inverse</code> as well.
template <class Preconditioner>
class SchurComplement : public Subscriptor
{
public:
SchurComplement (const BlockSparseMatrix<double> &A,
- const InverseMatrix<SparseMatrix<double>, Preconditioner> &Minv);
+ const InverseMatrix<SparseMatrix<double>, Preconditioner> &Minv);
void vmult (Vector<double> &dst,
- const Vector<double> &src) const;
+ const Vector<double> &src) const;
private:
const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
template <class Preconditioner>
SchurComplement<Preconditioner>::
SchurComplement (const BlockSparseMatrix<double> &A,
- const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv)
- :
- system_matrix (&A),
- m_inverse (&Minv),
- tmp1 (A.block(0,0).m()),
- tmp2 (A.block(0,0).m())
+ const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv)
+ :
+ system_matrix (&A),
+ m_inverse (&Minv),
+ tmp1 (A.block(0,0).m()),
+ tmp2 (A.block(0,0).m())
{}
template <class Preconditioner>
void SchurComplement<Preconditioner>::vmult (Vector<double> &dst,
- const Vector<double> &src) const
+ const Vector<double> &src) const
{
system_matrix->block(0,1).vmult (tmp1, src);
m_inverse->vmult (tmp2, tmp1);
}
- // @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 <code>degree+1</code> for
- // the vector-valued velocity components and
- // of order <code>degree</code> 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
- // <code>Triangulation::MeshSmoothing</code>
- // 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
+ // <code>degree+1</code> for the
+ // vector-valued velocity components and of
+ // order <code>degree</code> 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
+ // <code>Triangulation::MeshSmoothing</code>
+ // for details.
template <int dim>
StokesProblem<dim>::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 <int dim>
void StokesProblem<dim>::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
- // <code>DoFRenumbering::Cuthill_McKee</code>,
- // and then we renumber once again by
- // components. Since
- // <code>DoFRenumbering::component_wise</code>
- // 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
+ // <code>DoFRenumbering::Cuthill_McKee</code>,
+ // and then we renumber once again by
+ // components. Since
+ // <code>DoFRenumbering::component_wise</code>
+ // 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 <code>dim</code>
- // blocks for the velocity components,
- // so they can all be grouped in only
- // one block. The vector
- // <code>block_component</code> 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 <code>dim</code>
+ // blocks for the velocity components, so
+ // they can all be grouped in only one
+ // block. The vector
+ // <code>block_component</code> does
+ // precisely this: velocity values
+ // correspond to block 0, and pressure
+ // values will sit in block 1.
std::vector<unsigned int> 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
- // <code>DoFTools::count_dofs_per_block</code>
- // does the same as
- // <code>DoFTools::count_dofs_per_component</code>,
- // but now grouped as velocity and
- // pressure block via <code>block_component</code>.
+ // 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
+ // <code>DoFTools::count_dofs_per_block</code>
+ // does the same as
+ // <code>DoFTools::count_dofs_per_component</code>,
+ // but now grouped as velocity and pressure
+ // block via <code>block_component</code>.
std::vector<unsigned int> 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_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
- // <code>DoFTools::max_couplings_between_dofs</code>
- // 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 <code>csp</code> will be
- // released once the information has been
- // copied to
- // <code>sparsity_pattern</code>.
+ // 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
+ // <code>DoFTools::max_couplings_between_dofs</code>
+ // 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
+ // <code>csp</code> will be released once
+ // the information has been copied to
+ // <code>sparsity_pattern</code>.
{
BlockCompressedSparsityPattern csp;
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);
}
- // @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 <int dim>
void StokesProblem<dim>::assemble_system ()
{
std::vector<unsigned int> 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<dim> pressure_boundary_values;
std::vector<double> boundary_values (n_face_q_points);
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
- // This starts the loop over all
- // cells. With the abbreviations
- // <code>extract_u</code> etc.
- // introduced above, it is
- // evident what is going on.
+ // This starts the loop over all
+ // cells. With the abbreviations
+ // <code>extract_u</code> etc.
+ // introduced above, it is
+ // evident what is going on.
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
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<dofs_per_cell; ++j)
- {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
const SymmetricTensor<2,dim>
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
- // <code> phi_i_p * phi_j_p </code>,
- // 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
- // <code>phi_i_p * phi_j_p</code>
- // is only non-zero when all the
- // other terms vanish (and the other
- // way around).
+ // Note the way we write
+ // the contributions <code>
+ // phi_i_p * phi_j_p
+ // </code>, 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
+ // <code>phi_i_p *
+ // phi_j_p</code> is only
+ // non-zero when all the
+ // other terms vanish (and
+ // the other way around).
//
// Note also that operator*
// is overloaded for
// 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<GeometryInfo<dim>::faces_per_cell;
++face_no)
}
}
- // 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<dofs_per_cell; ++i)
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
- // <code>component_mask</code> 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 <code>component_mask</code> 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);
system_rhs);
}
- // Before we're going to solve
- // this linear system, we generate
- // a preconditioner for the
- // velocity-velocity matrix,
- // i.e., <code>block(0,0)</code>
- // in the system matrix. As mentioned
- // above, this depends on the
- // spatial dimension. Since this
- // handled automatically by the
- // template <code><dim></code>
- // in <code>InnerPreconditioner</code>,
- // 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.,
+ // <code>block(0,0)</code> in the system
+ // matrix. As mentioned above, this depends
+ // on the spatial dimension. Since this
+ // handled automatically by the template
+ // <code><dim></code> in
+ // <code>InnerPreconditioner</code>, 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<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
A_preconditioner->initialize (system_matrix.block(0,0),
- typename InnerPreconditioner<dim>::type::AdditionalData());
+ typename InnerPreconditioner<dim>::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
- // <code>solve</code> function is
- // rather straigt-forward and done in
- // a similar way as in step-20. To
- // start with, we need an object of
- // the <code>InverseMatrix</code> 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 <code>InnerPreconditioner<dim></code>.
+ // After the discussion in the introduction
+ // and the definition of the respective
+ // classes above, the implementation of the
+ // <code>solve</code> function is rather
+ // straigt-forward and done in a similar way
+ // as in step-20. To start with, we need an
+ // object of the <code>InverseMatrix</code>
+ // 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
+ // <code>InnerPreconditioner<dim></code>.
template <int dim>
void StokesProblem<dim>::solve ()
{
A_inverse (system_matrix.block(0,0), *A_preconditioner);
Vector<double> 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<double> schur_rhs (solution.block(1).size());
A_inverse.vmult (tmp, system_rhs.block(0));
SchurComplement<typename InnerPreconditioner<dim>::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
- // <code>InverseMatrix</code>, which
- // implements the <code>vmult</code>
- // 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
+ // <code>InverseMatrix</code>, which
+ // implements the <code>vmult</code>
+ // 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<SparseMatrix<double>,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;
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 -
- // <code>velocity</code> to the dim
- // components of velocity and
- // <code>p</code> 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
- // <code>DataComponentInterpretation</code>
- // 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
+ // - <code>velocity</code> to the dim
+ // components of velocity and <code>p</code>
+ // 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
+ // <code>DataComponentInterpretation</code>
+ // class. The rest of the function is then
+ // the same as in step-20.
template <int dim>
void
StokesProblem<dim>::output_results (const unsigned int refinement_cycle) const
= DataComponentInterpretation::component_is_part_of_vector;
data_out.add_data_vector (solution, solution_names,
- DataOut<dim>::type_dof_data,
- data_component_interpretation);
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
data_out.build_patches ();
}
- // @sect4{StokesProblem::refine_mesh}
+ // @sect4{StokesProblem::refine_mesh}
- // This is the last interesting function
- // of the <code>StokesProblem</code> 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 <code>StokesProblem</code> 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 <int dim>
void
StokesProblem<dim>::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 <int dim>
void StokesProblem<dim>::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<unsigned int> subdivisions (dim, 1);
subdivisions[0] = 4;
Point<dim>(2,0) :
Point<dim>(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<dim>::active_cell_iterator
- cell = triangulation.begin_active();
+ cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
for (unsigned int f=0; f<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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)
{
}
- // @sect3{The <code>main</code> function}
+ // @sect3{The <code>main</code> 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
projects / dealii-svn.git / commitdiff