/* $Id$ */
/* */
-/* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors */
+/* Copyright (C) 2010 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
// The first few files have already
// been covered in step-12
// and will thus not be further
- // commented on.
+ // commented on:
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <lac/vector.h>
#include <base/timer.h>
// Here come the new include files
- // for using the MeshWorker framework
+ // for using the MeshWorker framework:
#include <numerics/mesh_worker.h>
#include <numerics/mesh_worker_loop.h>
// @sect3{Equation data}
//
+ // First, we need to describe the
+ // coefficients in the equation. Here, this
+ // concerns the boundary values, which we
+ // choose in the same way as for step-12:
template <int dim>
class BoundaryValues: public Function<dim>
{
const unsigned int component=0) const;
};
- // The inflow boundary of the
- // unit square [0,1]^2 are the right
- // and the lower boundaries. We
- // prescribe discontinuous boundary
- // values 1 and 0 on the x-axis and
- // value 0 on the right boundary. The
- // values of this function on the
- // outflow boundaries will not be
- // used within the DG scheme.
+ // Given the flow direction, the inflow
+ // boundary of the unit square $[0,1]^2$ are
+ // the right and the lower boundaries. We
+ // prescribe discontinuous boundary values 1
+ // and 0 on the x-axis and value 0 on the
+ // right boundary. The values of this
+ // function on the outflow boundaries will
+ // not be used within the DG scheme.
template <int dim>
void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
std::vector<double> &values,
// the MeshWorker framework. Since it
// will be used by
// MeshWorker::AssemblingIntegrator,
- // it needs the functions for cell,
+ // it needs functions for cell,
// boundary and interior face
// integration specified exactly as
// below.
beta(1) = fe_v.quadrature_point(point)(0);
beta /= beta.norm();
+ // We solve a homogeneous
+ // equation, thus no right
+ // hand side shows up in
+ // the cell term.
+ // What's left is
+ // integrating the matrix entries.
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- {
- // We solve a homogeneous
- // equation, thus no right
- // hand side shows up in
- // the cell term.
- // What's left is
- // integrating the matrix entries.
- for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
- local_matrix(i,j) -= beta*fe_v.shape_grad(i,point)*
- fe_v.shape_value(j,point) *
- JxW[point];
- }
+ for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+ local_matrix(i,j) -= beta*fe_v.shape_grad(i,point)*
+ fe_v.shape_value(j,point) *
+ JxW[point];
}
}
- // Now the same for the boundary
- // terms. Note that now we use
- // FEFaceValuesBase in order to get
+ // Now the same for the boundary terms. Note
+ // that now we use FEFaceValuesBase, the base
+ // class for both FEFaceValues and
+ // FESubfaceValues, in order to get access to
// normal vectors.
template <int dim>
void DGIntegrator<dim>::bdry(FaceInfo& info) const
// four local matrices. The letters
// u and v refer to trial and test
// functions, respectively. The
- // numbers indicate the cells
+ // %numbers indicate the cells
// provided by info1 and info2. By
// convention, the two matrices in
// each info object refer to the
// vectors. Fortunately, the
// interface terms only involve the
// solution and the right hand side
- // does not obtain a contribution.
+ // does not receive any contributions.
const std::vector<double> &JxW = fe_v.get_JxW_values ();
const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
if (beta_n>0)
{
// This term we've already
- // seen.
+ // seen:
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
u1_v1_matrix(i,j) += beta_n *
else
{
// This one we've already
- // seen, too.
+ // seen, too:
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
u2_v1_matrix(i,l) += beta_n *
// And this is another new
// one: $(\beta\cdot n \hat
// u,\hat v)_{\partial
- // \kappa_-}$.
+ // \kappa_-}$:
for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
u2_v2_matrix(k,l) -= beta_n *
// use a DG method of a different
// degree the whole program stays
// the same, only replace 1 in
- // the constructor by the wanted
- // degree.
+ // the constructor by the desired
+ // polynomial degree.
FE_DGQ<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
- // And there are two solution
- // vectors, that store the
- // solutions to the problems
- // corresponding to the two
+ // In step-12 we had two solution vectors
+ // that stored the solutions to the
+ // problems corresponding to the two
// different assembling routines
// <code>assemble_system1</code> and
- // <code>assemble_system2</code>;
+ // <code>assemble_system2</code>. In this
+ // program, the goal is only to show the
+ // MeshWorker framework, so we only
+ // assemble the system in one of the two
+ // ways, and consequently we have only
+ // one solution vector along with the
+ // single <code>assemble_system</code>
+ // function declared above:
Vector<double> solution;
Vector<double> right_hand_side;
};
+ // We start with the
+ // constructor. This is the
+ // place to change the
+ // polynomial degree of the
+ // finite element shape
+ // functions.
template <int dim>
DGMethod<dim>::DGMethod ()
:
- // Change here for DG
- // methods of
- // different degrees.
fe (1),
dof_handler (triangulation)
{}
}
+ // In the function that sets up the usual
+ // finite element data structures, we first
+ // need to distribute the DoFs.
template <int dim>
void DGMethod<dim>::setup_system ()
{
- // First we need to distribute the
- // DoFs.
dof_handler.distribute_dofs (fe);
- // The DoFs of a cell are coupled
- // with all DoFs of all neighboring
- // cells. Therefore the maximum
- // number of matrix entries per row
- // is needed when all neighbors of
- // a cell are once more refined
- // than the cell under
+ // The DoFs of a cell are coupled with all
+ // DoFs of all neighboring cells, along
+ // with all of its siblings on the current
+ // cell. Therefore the maximum number of
+ // matrix entries per row is needed when
+ // all neighbors of a cell are once more
+ // refined than the cell under
// consideration.
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
- (GeometryInfo<dim>::faces_per_cell
- *GeometryInfo<dim>::max_children_per_face+1)*fe.dofs_per_cell);
+ (GeometryInfo<dim>::faces_per_cell *
+ GeometryInfo<dim>::max_children_per_face
+ +
+ 1)*fe.dofs_per_cell);
- // For DG discretizations we call
- // the function analogue to
- // DoFTools::make_sparsity_pattern.
+ // To build the sparsity pattern for DG
+ // discretizations, we can call the
+ // function analogue to
+ // DoFTools::make_sparsity_pattern, which
+ // is called
+ // DoFTools::make_flux_sparsity_pattern:
DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
// All following function calls are
}
// @sect4{Function: assemble_system}
- // Here we see the major difference
- // to assembling by hand. Instead of
- // writing loops over cells and
- // faces, we leave all this to the
- // MeshWorker framework. In order to
- // do so, we just have to define
- // local integration objects and use
- // one of the classes in Assembler to
- // build the global system.
+
+ // Here we see the major difference to
+ // assembling by hand. Instead of writing
+ // loops over cells and faces, we leave all
+ // this to the MeshWorker framework. In order
+ // to do so, we just have to define local
+ // integration objects and use one of the
+ // classes in namespace MeshWorker::Assembler
+ // to build the global system.
template <int dim>
void DGMethod<dim>::assemble_system ()
{
// the global sparse matrix and the
// right hand side vector.
//
- // It turns out,
// MeshWorker::AssemblingIntegrator
- // itself is not that clever at
- // all, but all its capabilities
- // are provided by the two later
+ // is not all that clever by itself,
+ // but its capabilities
+ // are provided by its two latter
// template arguments. By
// exchanging
// MeshWorker::Assembler::SystemSimple,
// we could for instance assemble a
// BlockMatrix or just a Vector
// instead.
- MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> >, DGIntegrator<dim> >
+ MeshWorker::AssemblingIntegrator
+ <dim,
+ MeshWorker::Assembler::SystemSimple<SparseMatrix<double>,
+ Vector<double> >,
+ DGIntegrator<dim> >
integrator(dg);
// First, we initialize the
// independently, we have to hand
// over this value three times.
const unsigned int n_gauss_points = dof_handler.get_fe().degree+1;
- integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
+ integrator.initialize_gauss_quadrature(n_gauss_points,
+ n_gauss_points,
+ n_gauss_points);
// These are the types of values we
// need for integrating our
// and interior faces, as well as
// interior neighbor faces, which is
// forced by the four @p true values.
- UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients;
+ UpdateFlags update_flags = update_quadrature_points |
+ update_values |
+ update_gradients;
integrator.add_update_flags(update_flags, true, true, true, true);
// Finally, we have to tell the
- // assembler base class, where to
+ // assembler base class where to
// put the local data. These will
// be our system matrix and the
// right hand side.
integrator.initialize(system_matrix, right_hand_side);
- // Finally, we get to the
+ // We are now ready to get to the
// integration loop. @p info_box is
- // an object, that generates the
+ // an object that generates the
// extended iterators for cells and
// faces of type
// MeshWorker::IntegrationInfo. Since
//
// For this simple problem we use the
// simplest possible solver, called
- // Richardson iteration, that
- // represents a simple defect
- // correction. This, in combination
- // with a block SSOR preconditioner,
- // that uses the special block matrix
- // structure of system matrices
- // arising from DG
- // discretizations. The size of these
- // blocks are the number of DoFs per
- // cell. Here, we use a SSOR
- // preconditioning as we have not
- // renumbered the DoFs according to
- // the flow field. If the DoFs are
- // renumbered downstream the flow,
- // then a block Gauss-Seidel
+ // Richardson iteration, that represents a
+ // simple defect correction. This, in
+ // combination with a block SSOR
+ // preconditioner, that uses the special
+ // block matrix structure of system matrices
+ // arising from DG discretizations. The size
+ // of these blocks are the number of DoFs per
+ // cell. Here, we use a SSOR preconditioning
+ // as we have not renumbered the DoFs
+ // according to the flow field. If the DoFs
+ // are renumbered in the downstream direction
+ // of the flow, then a block Gauss-Seidel
// preconditioner (see the
// PreconditionBlockSOR class with
- // relaxation=1) makes a much better
- // job.
+ // relaxation=1) does a much better job.
template <int dim>
void DGMethod<dim>::solve (Vector<double> &solution)
{
// preconditioner,
PreconditionBlockSSOR<SparseMatrix<double> > preconditioner;
- // we assigned the matrix to it and
- // set the right block size.
+ // then assign the matrix to it and
+ // set the right block size:
preconditioner.initialize(system_matrix, fe.dofs_per_cell);
// After these preparations we are
// The following <code>run</code> function is
- // similar to previous examples. The
- // only difference is that the
- // problem is assembled and solved
- // twice on each refinement step;
- // first by <code>assemble_system1</code> that
- // implements the first version and
- // then by <code>assemble_system2</code> that
- // implements the second version of
- // writing the DG
- // discretization. Furthermore the
- // time needed by each of the two
- // assembling routines is measured.
+ // similar to previous examples.
template <int dim>
void DGMethod<dim>::run ()
{
assemble_system ();
solve (solution);
- // Finally we perform the
- // output.
output_results (cycle);
}
}
// The following <code>main</code> function is
- // similar to previous examples and
- // need not to be commented on.
+ // similar to previous examples as well, and
+ // need not be commented on.
int main ()
{
try