-/* $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 */
/* further information on this license. */
-
+ // @sect3{Include 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 <lac/block_sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
-#include <lac/sparse_direct.h>
-#include <lac/sparse_ilu.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <dofs/dof_constraints.h>
#include <fe/fe_q.h>
-#include <fe/fe_dgq.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <fe/mapping_q1.h>
-#include <fe/mapping_c1.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
-#include <numerics/solution_transfer.h>
+ // 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.
+#include <lac/sparse_ilu.h>
+
+ // This is C++:
#include <fstream>
#include <sstream>
+ // 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 <int dim>
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 <code>SparseILU@<double></code>.
template <>
struct InnerPreconditioner<3>
{
};
-
+ // @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.
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.
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>.
+
+ // 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>
{
+ // We implement similar functions
+ // for the right hand side.
template <int dim>
class RightHandSide : public Function<dim>
{
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
+
};
+ // @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
+ // <code>extract_grad_s_u</code>, 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 <int dim>
Tensor<1,dim>
extract_u (const FEValuesBase<dim> &fe_values,
}
-
template <int dim>
Tensor<2,dim>
extract_grad_s_u (const FEValuesBase<dim> &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 <int dim>
+double
+scalar_product (const Tensor<2,dim> &a,
+ const Tensor<2,dim> &b)
+{
+ double tmp = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ tmp += a[i][j] * b[i][j];
+ return tmp;
+}
-
+
+
+ // @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.
+
+ // @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.
template <class Matrix, class Preconditioner>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix (const Matrix &m,
- const Preconditioner &preconditioner);
+ const Preconditioner &preconditioner);
void vmult (Vector<double> &dst,
- const Vector<double> &src) const;
+ const Vector<double> &src) const;
private:
const SmartPointer<const Matrix> matrix;
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:
+
+ // 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.
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.
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}
+
+ // @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 discussion of
+ // <code>Triangulation::MeshSmoothing</code>
+ // for details.
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int degree)
:
degree (degree),
- triangulation (Triangulation<dim>::maximum_smoothing),
+ triangulation (Triangulation<dim>::maximum_smoothing),
fe (FE_Q<dim>(degree+1), dim,
FE_Q<dim>(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 <int dim>
void StokesProblem<dim>::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
+ // <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);
- DoFRenumbering::component_wise (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.
+ 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.
hanging_node_constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
- hanging_node_constraints);
+ hanging_node_constraints);
hanging_node_constraints.close ();
- std::vector<unsigned int> 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
+ // <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_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
+ // <code>max_couplings_between_dofs@</code>
+ // yields a very large number for the
+ // coupling between the individual dofs,
+ // so that the memory initially provided
+ // in the <code>reinit</code> 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
+ // <code>BlockCompressedSparsityPattern</code>,
+ // 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;
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);
}
-template <int dim>
-double
-scalar_product (const Tensor<2,dim> &a,
- const Tensor<2,dim> &b)
-{
- double tmp = 0;
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- tmp += a[i][j] * b[i][j];
- return tmp;
-}
-
-
+ // @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.
template <int dim>
void StokesProblem<dim>::assemble_system ()
QGauss<dim-1> face_quadrature_formula(degree+2);
FEValues<dim> 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<dim> 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;
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.
const PressureBoundaryValues<dim> pressure_boundary_values;
std::vector<double> boundary_values (n_face_q_points);
-
+
+ const RightHandSide<dim> right_hand_side;
+ std::vector<Vector<double> > rhs_values (n_q_points,
+ Vector<double>(dim+1));
+
+ // This starts the loop over all
+ // cells. With the abbreviations
+ // <code>extract_u</code> etc.
+ // introduced above, it is very
+ // clear what is going on.
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
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<n_q_points; ++q)
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
-
- const Tensor<1,dim> phi_i_u = extract_u (fe_values, i, q);
-
+ {
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
const Tensor<2,dim> 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<dofs_per_cell; ++j)
+ {
+ const Tensor<2,dim> 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<dofs_per_cell; ++j)
- {
- const Tensor<2,dim> 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
+ // <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.
+ 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<GeometryInfo<dim>::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);
.value_list (fe_face_values.get_quadrature_points(),
boundary_values);
- for (unsigned int q=0; q<n_face_q_points; ++q)
+ for (unsigned int q=0; q<n_face_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- 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));
- }
+ {
+ 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; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- local_matrix(i,j));
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
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
+ // componenent, so that the condensation
+ // is performed only on
+ // velocity dofs.
hanging_node_constraints.condense (system_matrix);
hanging_node_constraints.condense (system_rhs);
std::vector<bool> component_mask (dim+1, true);
component_mask[dim] = false;
VectorTools::interpolate_boundary_values (dof_handler,
- 1,
- BoundaryValues<dim>(),
- boundary_values,
- component_mask);
+ 1,
+ BoundaryValues<dim>(),
+ 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., <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());
- 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
+ // <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 ()
{
- const InverseMatrix<SparseMatrix<double>,typename InnerPreconditioner<dim>::type>
+ const InverseMatrix<SparseMatrix<double>,
+ typename InnerPreconditioner<dim>::type>
A_inverse (system_matrix.block(0,0), *A_preconditioner);
Vector<double> tmp (solution.block(0).size());
- Vector<double> 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<double> 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<typename InnerPreconditioner<dim>::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
+ // <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
+ // 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<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);
+ 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 -
+ // <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
std::vector<std::string> solution_names (dim, "velocity");
solution_names.push_back ("p");
- DataOut<dim> data_out;
+ DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
= 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 ();
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 <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 ()
std::vector<bool> component_mask (dim+1, false);
component_mask[dim] = true;
KellyErrorEstimator<dim>::estimate (dof_handler,
- QGauss<dim-1>(degree+1),
- typename FunctionMap<dim>::type(),
- solution,
- estimated_error_per_cell,
- component_mask);
+ QGauss<dim-1>(degree+1),
+ typename FunctionMap<dim>::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 <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.
std::vector<unsigned int> subdivisions (dim, 1);
subdivisions[0] = 4;
-
+
GridGenerator::subdivided_hyper_rectangle (triangulation,
- subdivisions,
- (dim == 2 ?
- Point<dim>(-2,-1) :
- Point<dim>(-2,0,-1)),
- (dim == 2 ?
- Point<dim>(2,0) :
- Point<dim>(2,1,0)));
+ subdivisions,
+ (dim == 2 ?
+ Point<dim>(-2,-1) :
+ Point<dim>(-2,0,-1)),
+ (dim == 2 ?
+ 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.
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.
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 <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.
int main ()
{
try
projects / dealii-svn.git / commitdiff