/* $Id$ */
-/* Author: Ralf Hartmann, University of Heidelberg, 2000 */
+/* Author: Ralf Hartmann, University of Heidelberg, 2001 */
// The first few files have already
// been covered in previous examples
// used as all other finite elements.
#include <fe/fe_dgq.h>
// We are going to use the simplest
- // possible solver, called richardson
+ // possible solver, called Richardson
// iteration, that represents a simple
// defect correction. This, in
// combination with a block SSOR
#include <fstream>
- // First we define the class
+ // @sect3{Equation data}
+ //
+ // First we define the classes
// representing the equation-specific
// functions. Both classes, ``RHS''
// and ``BoundaryValues'', are
- // derived from the Function
+ // derived from the ``Function''
// class. Only the ``value_list''
// function are implemented because
// only lists of function values are
class RHS: public Function<dim>
{
public:
- RHS() {};
-
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component=0) const;
class BoundaryValues: public Function<dim>
{
public:
- BoundaryValues() {};
-
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component=0) const;
};
- // The class ``Beta'' that represents
- // the vector valued flow field of
- // the linear transport equation is
- // not derived from the Function
+ // The class ``Beta'' represents the
+ // vector valued flow field of the
+ // linear transport equation and is
+ // not derived from the ``Function''
// class as we prefer to get function
// values of type ``Point'' rather
// than of type
class Beta
{
public:
- Beta () {};
-
void value_list (const std::vector<Point<dim> > &points,
std::vector<Point<dim> > &values) const;
};
// ``value_list'' functions of these
// classes are rather simple. For
// simplicity the right hand side is
- // set to be zero.
+ // set to be zero but will be
+ // assembled anyway.
template <int dim>
-void RHS<dim>::value_list(const std::vector<Point<dim> > &,
+void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int) const
{
+ // Usually we check whether input
+ // parameter have the right sizes.
+ Assert(values.size()==points.size(),
+ ExcDimensionMismatch(values.size(),points.size()));
+
for (unsigned int i=0; i<values.size(); ++i)
values[i]=0;
}
// The flow field is chosen to be
- // circular, anticlockwise, and with
+ // circular, counterclockwise, and with
// the origin as midpoint.
-template <>
-void Beta<2>::value_list(const std::vector<Point<2> > &points,
- std::vector<Point<2> > &values) const
+template <int dim>
+void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
+ std::vector<Point<dim> > &values) const
{
Assert(values.size()==points.size(),
ExcDimensionMismatch(values.size(),points.size()));
for (unsigned int i=0; i<points.size(); ++i)
{
- const Point<2> &p=points[i];
- Point<2> &beta=values[i];
+ const Point<dim> &p=points[i];
+ Point<dim> &beta=values[i];
beta(0)=-p(1);
beta(1)=p(0);
- beta/=sqrt(beta*beta);
+ beta/=sqrt(beta.square());
}
}
// @sect3{Class: DGTransportEquation}
//
- // Next we define the equation-
- // dependent and DG-method-dependent
- // class ``DGTransportEquation''. Its
+ // Next we define the
+ // equation-dependent and
+ // DG-method-dependent class
+ // ``DGTransportEquation''. Its
// member functions were already
// mentioned in the Introduction and
// will be explained
void assemble_cell_term(const FEValues<dim>& fe_v,
FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector);
+ Vector<double> &cell_vector) const;
void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector);
+ Vector<double> &cell_vector) const;
void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
const FEFaceValuesBase<dim>& fe_v_neighbor,
FullMatrix<double> &u_v_matrix,
- FullMatrix<double> &un_v_matrix);
+ FullMatrix<double> &un_v_matrix) const;
void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
const FEFaceValuesBase<dim>& fe_v_neighbor,
FullMatrix<double> &u_v_matrix,
FullMatrix<double> &un_v_matrix,
FullMatrix<double> &u_vn_matrix,
- FullMatrix<double> &un_vn_matrix);
+ FullMatrix<double> &un_vn_matrix) const;
private:
- Beta<dim> beta_function;
- RHS<dim> rhs_function;
- BoundaryValues<dim> boundary_function;
+ const Beta<dim> beta_function;
+ const RHS<dim> rhs_function;
+ const BoundaryValues<dim> boundary_function;
};
// it is of size 4 times 4, and
// ``cell_vector'' is of size 4.
// When this function is invoked,
- // ``fe_v'' was reinited with the
+ // ``fe_v'' is already reinit'ed with the
// current cell before and includes
// all shape values needed.
template <int dim>
void DGTransportEquation<dim>::assemble_cell_term(
- const FEValues<dim>& fe_v,
+ const FEValues<dim> &fe_v,
FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector)
+ Vector<double> &cell_vector) const
{
// First we ask ``fe_v'' for the
// shape gradients, shape values and
// quadrature weights,
- const vector<vector<Tensor<1,2> > > &grad_v = fe_v.get_shape_grads ();
+ const std::vector<std::vector<Tensor<1,dim> > > &grad_v = fe_v.get_shape_grads ();
const FullMatrix<double> &v = fe_v.get_shape_values ();
- const vector<double> &JxW = fe_v.get_JxW_values ();
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();
// Then the flow field beta and the
// ``rhs_function'' are evaluated at
// the quadrature points,
- vector<Point<dim> > beta (fe_v.n_quadrature_points);
- vector<double> rhs (fe_v.n_quadrature_points);
+
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
+ std::vector<double> rhs (fe_v.n_quadrature_points);
beta_function.value_list (fe_v.get_quadrature_points(), beta);
rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
// function assembles the face terms
// at boundary faces. When this
// function is invoked, ``fe_v'' is
- // already reinited with the current
+ // already reinit'ed with the current
// cell and current face. Hence it
// provides the shape values on that
// boundary face.
void DGTransportEquation<dim>::assemble_boundary_term(
const FEFaceValues<dim>& fe_v,
FullMatrix<double> &u_v_matrix,
- Vector<double> &cell_vector)
+ Vector<double> &cell_vector) const
{
// First we check whether the
// current face is really at the
Assert(fe_v.get_face()->at_boundary(), ExcInternalError());
// Again, as in the previous
- // function, we ask the FEValues
+ // function, we ask the ``FEValues''
// object for the shape values and
// the quadrature weights
const FullMatrix<double> &v = fe_v.get_shape_values ();
- const vector<double> &JxW = fe_v.get_JxW_values ();
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();
// but here also for the normals.
- const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const
std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
// We evaluate the flow field
// and the boundary values at the
// quadrature points.
- vector<Point<dim> > beta (fe_v.n_quadrature_points);
- vector<double> g(fe_v.n_quadrature_points);
+
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
+ std::vector<double> g(fe_v.n_quadrature_points);
beta_function.value_list (fe_v.get_quadrature_points(), beta);
boundary_function.value_list (fe_v.get_quadrature_points(), g);
// introduction.
for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
{
- double beta_n=beta[point] * normals[point];
- if (beta_n>0)
+ const double beta_n=beta[point] * normals[point];
// We assemble the term
// $(\beta\cdot n
// u,v)_{\partial K_+}$,
+ if (beta_n>0)
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
u_v_matrix(i,j) += beta_n *
//
// When this function is invoked,
// ``fe_v'' and ``fe_v_neighbor'' are
- // already reinited with the current
+ // already reinit'ed with the current
// cell and the neighoring cell,
// respectively, as well as with the
// current face. Hence they provide
// contributions to the system matrix
// that are based on outer values of
// u, see $\hat u_h$ in the
- // Introduction, and inner values of
+ // introduction, and inner values of
// v, see $v_h$. Here we note that
// ``un'' is the short notation for
// ``u_neighbor'' and represents
const FEFaceValuesBase<dim>& fe_v,
const FEFaceValuesBase<dim>& fe_v_neighbor,
FullMatrix<double> &u_v_matrix,
- FullMatrix<double> &un_v_matrix)
+ FullMatrix<double> &un_v_matrix) const
{
// First we check that the current
// face is not at the boundary by
// normals
const FullMatrix<double> &v = fe_v.get_shape_values ();
const FullMatrix<double> &v_neighbor = fe_v_neighbor.get_shape_values ();
- const vector<double> &JxW = fe_v.get_JxW_values ();
- const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();
+ const
std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
// and we evaluate the flow field
// at the quadrature points.
- vector<Point<dim> > beta (fe_v.n_quadrature_points);
+
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
beta_function.value_list (fe_v.get_quadrature_points(), beta);
// Then we assemble the cell
// introduction.
for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
{
- double beta_n=beta[point] * normals[point];
- if (beta_n>0)
+ const double beta_n=beta[point] * normals[point];
// We assemble the term
// $(\beta\cdot n
// u,v)_{\partial K_+}$,
+ if (beta_n>0)
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
u_v_matrix(i,j) += beta_n *
FullMatrix<double> &u_v_matrix,
FullMatrix<double> &un_v_matrix,
FullMatrix<double> &u_vn_matrix,
- FullMatrix<double> &un_vn_matrix)
+ FullMatrix<double> &un_vn_matrix) const
{
// the first few lines are the same
const FullMatrix<double> &v = fe_v.get_shape_values ();
const FullMatrix<double> &v_neighbor = fe_v_neighbor.get_shape_values ();
- const vector<double> &JxW = fe_v.get_JxW_values ();
- const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();
+ const
std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
- vector<Point<dim> > beta (fe_v.n_quadrature_points);
+
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
beta_function.value_list (fe_v.get_quadrature_points(), beta);
for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
{
- double beta_n=beta[point] * normals[point];
+ const double beta_n=beta[point] * normals[point];
if (beta_n>0)
{
// This terms we've already seen.
// the main class of step 6. One of
// the differences is that there's no
// ConstraintMatrix object. This is,
- // because there are no hanging nodes
- // in DG discretizations.
+ // because there are no hanging node
+ // constraints in DG discretizations.
template <int dim>
class DGMethod
{
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
- MappingQ1<dim> mapping;
+ const MappingQ1<dim> mapping;
- // Furthermore we want to
- // use DG elements of degree 1
- // (but this is only specified in
- // the constructor):
+ // Furthermore we want to use DG
+ // elements of degree 1 (but this
+ // is only specified in the
+ // constructor). If you want to
+ // use a DG method of a different
+ // degree the whole program stays
+ // the same, only replace 1 in
+ // the constructor by the wanted
+ // degree.
FE_DGQ<dim> fe;
DoFHandler<dim> dof_handler;
// formulae for the cell and the
// face terms of the
// discretization.
- QGauss4<dim> quadrature;
- QGauss4<dim-1> face_quadrature;
+ const QGauss4<dim> quadrature;
+ const QGauss4<dim-1> face_quadrature;
// And there are two solution
// vectors, that store the
// object of the
// DGTransportEquations class
// described above.
- DGTransportEquation<dim> dg;
+ const DGTransportEquation<dim> dg;
};
- // Now for the implementation of the
- // main class. Constructor and
- // destructor follow the same
- // pattern that was used previously,
- // so we need not comment on these
- // two functions:
template <int dim>
DGMethod<dim>::DGMethod () :
+ // Change here for DG
+ // methods of
+ // different degrees.
fe (1),
dof_handler (triangulation)
{}
// The DoFs of a cell are coupled
// with all DoFs of all neighboring
// cells. Therefore the maximum
- // number of matrix entries is
- // needed when all neighbors of a
- // cell are once more refined than
- // the cell under consideration.
+ // 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
// ``assemble_cell_term'',
// ``assemble_boundary_term'' and
// ``assemble_face_term1'' functions
- // of the DGTransportEquation object.
- // The ``assemble_face_term1''
- // function takes two
- // FEFaceValuesBase objects; one for
- // the shape functions on the current
- // cell and the other for shape
- // functions on the neighboring cell
- // under consideration. Both objects
- // are either of class FEFaceValues
- // or of class FESubfaceValues (both
- // derived from FEFaceValuesBase)
- // according to following cases
- // already mentioned in the
- // introduction:
+ // of the ``DGTransportEquation''
+ // object. The
+ // ``assemble_boundary_term'' covers
+ // the first case mentioned in the
+ // introduction.
+ //
+ // 1. face is at boundary
//
- // 1. face is at boundary (current
- // cell: FEFaceValues);
+ // This function takes a
+ // ``FEFaceValues'' object as
+ // argument. In contrast to that
+ // ``assemble_face_term1''
+ // takes two ``FEFaceValuesBase''
+ // objects; one for the shape
+ // functions on the current cell and
+ // the other for shape functions on
+ // the neighboring cell under
+ // consideration. Both objects are
+ // either of class ``FEFaceValues''
+ // or of class ``FESubfaceValues''
+ // (both derived from
+ // ``FEFaceValuesBase'') according to
+ // the remaining cases mentioned
+ // in the introduction:
//
// 2. neighboring cell is finer
- // (current cell: FESubfaceValues,
- // neighboring cell: FEFaceValues);
+ // (current cell: ``FESubfaceValues'',
+ // neighboring cell: ``FEFaceValues'');
//
// 3. neighboring cell is of the same
// refinement level (both, current
// and neighboring cell:
- // FEFaceValues);
+ // ``FEFaceValues'');
//
// 4. neighboring cell is coarser
- // (current cell: FEFaceValues,
+ // (current cell: ``FEFaceValues'',
// neighboring cell:
- // FESubfaceValues).
+ // ``FESubfaceValues'').
//
// If we considered globally refined
- // meshes then only cases 1 and 3
- // would occur. But as we consider
- // also locally refined meshes we
- // need to distinguish all four cases
- // making the following assembling
- // function a bit longish.
+ // meshes then only case 3 would
+ // occur. But as we consider also
+ // locally refined meshes we need to
+ // distinguish all four cases making
+ // the following assembling function
+ // a bit longish.
template <int dim>
void DGMethod<dim>::assemble_system1 ()
{
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
- vector<unsigned int> dofs (dofs_per_cell);
- vector<unsigned int> dofs_neighbor (dofs_per_cell);
-
- // First we create the Update flags
- // for the FEValues and the
- // FEFaceValues objects.
- UpdateFlags update_flags = UpdateFlags(update_values
- | update_gradients
- | update_q_points
- | update_JxW_values);
-
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+ // First we create the
+ // ``UpdateFlags'' for the
+ // ``FEValues'' and the
+ // ``FEFaceValues'' objects.
+ UpdateFlags update_flags = update_values
+ | update_gradients
+ | update_q_points
+ | update_JxW_values;
+
// Note, that on faces we do not
// need gradients but we need
// normal vectors.
- UpdateFlags face_update_flags = UpdateFlags(update_values
- | update_q_points
- | update_JxW_values
- | update_normal_vectors);
+ UpdateFlags face_update_flags = update_values
+ | update_q_points
+ | update_JxW_values
+ | update_normal_vectors;
// On the neighboring cell we only
// need the shape values. Given a
// the normal vectors are known to
// be the negative of the normal
// vectors of the current cell.
- UpdateFlags neighbor_face_update_flags = UpdateFlags(update_values);
+ UpdateFlags neighbor_face_update_flags = update_values;
- // Then we create the FEValues
+ // Then we create the ``FEValues''
// object. Note, that since version
- // 3.2.0 the constructor of this
- // class takes a Mapping object as
- // first argument. Although the
- // constructor without Mapping
+ // 3.2.0 of deal.II the constructor
+ // of this class takes a
+ // ``Mapping'' object as first
+ // argument. Although the
+ // constructor without ``Mapping''
// argument is still supported it
// is recommended to use the new
// constructor. This reduces the
// effect of `hidden magic' (the
// old constructor implicitely
- // assumes a MappingQ1 mapping) and
- // makes it easier to change the
- // Mapping object later.
+ // assumes a ``MappingQ1'' mapping)
+ // and makes it easier to change
+ // the mapping object later.
FEValues<dim> fe_v (
mapping, fe, quadrature, update_flags);
// Similarly we create the
- // FEFaceValues and FESubfaceValues
- // objects for both, the current
- // and the neighboring cell. Within
- // the following nested loop over
- // all cells and all faces of the
- // cell they will be reinited to
- // the current cell and the face
- // (and subface) number.
+ // ``FEFaceValues'' and
+ // ``FESubfaceValues'' objects for
+ // both, the current and the
+ // neighboring cell. Within the
+ // following nested loop over all
+ // cells and all faces of the cell
+ // they will be reinited to the
+ // current cell and the face (and
+ // subface) number.
FEFaceValues<dim> fe_v_face (
mapping, fe, face_quadrature, face_update_flags);
FESubfaceValues<dim> fe_v_subface (
// include the `u and v terms' and
// the second that will include the
// `un and v terms'. Here we recall
- // our the convention that `un' is
- // the short cut for `u_neighbor'
+ // the convention that `un' is
+ // the shortcut for `u_neighbor'
// and represents the $u_hat$, see
// introduction.
FullMatrix<double> u_v_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_vector (dofs_per_cell);
// Furthermore we need some cell
- // and face iterators
- DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- DoFHandler<dim>::face_iterator face;
- DoFHandler<dim>::cell_iterator neighbor;
- DoFHandler<dim>::active_cell_iterator neighbor_child;
+ // iterators.
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
// Now we start the loop over all
- // active cells
+ // active cells.
for (;cell!=endc; ++cell)
{
- // and reinit the FEValues
- // object for the current cell,
+ // In the
+ // ``assemble_face_term1''
+ // function contributions to
+ // the cell matrices and the
+ // cell vector are only
+ // ADDED. Therefore on each
+ // cell we need to reset the
+ // ``u_v_matrix'' and
+ // ``cell_vector'' to zero,
+ // before assembling the cell terms.
+ u_v_matrix.clear ();
+ cell_vector.clear ();
+
+ // Now we reinit the ``FEValues''
+ // object for the current cell
fe_v.reinit (cell);
- // Call the function that
- // assembles the cell
+ // and call the function
+ // that assembles the cell
// terms. The first argument is
- // the FEValues that was
- // already reinited on current
- // the cell.
+ // the ``FEValues'' that was
+ // previously reinit'ed on the
+ // current cell.
dg.assemble_cell_term(fe_v,
u_v_matrix,
cell_vector);
- // As in previous example steps
- // the vector `dofs' includes
- // the dof_indices.
+ // As in previous examples the
+ // vector `dofs' includes the
+ // dof_indices.
cell->get_dof_indices (dofs);
// This is the start of the
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
// First we set the face
- // iterator.
- face = cell->face(face_no);
-
+ // iterator
+ typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+
+ // and clear the
+ // ``un_v_matrix'' on each
+ // face.
+ un_v_matrix.clear();
+
// Now we distinguish the
// four different cases in
// the ordering mentioned
if (face->at_boundary())
{
// We reinit the
- // FEFaceValues object
- // to the current face
+ // ``FEFaceValues''
+ // object to the
+ // current face
fe_v_face.reinit (cell, face_no);
// and assemble the
// domain, therefore
// there must exist a
// neighboring cell.
- neighbor = cell->neighbor(face_no);
+ typename DoFHandler<dim>::cell_iterator neighbor=
+ cell->neighbor(face_no);;
// We proceed with the
// second and most
// difference of not
// more than one, the
// neighboring cell is
- // known to be only
+ // known to be at most
// ONCE more refined
// than the current
// cell. Furthermore
// also the face is
- // once more refined,
+ // more refined,
// i.e. it has
- // children.
+ // children. Here we
+ // note that the
+ // following part of
+ // code will not work
+ // for ``dim==1''.
if (face->has_children())
{
- // first we store
+ // First we store
// which number the
// current cell has
// in the list of
// `behind' the
// current
// subface.
- neighbor_child = neighbor->child(GeometryInfo<dim>::
- child_cell_on_face(neighbor2,subface_no));
-
+ typename DoFHandler<dim>::active_cell_iterator neighbor_child=
+ neighbor->child(GeometryInfo<dim>::
+ child_cell_on_face(neighbor2,subface_no));
+
// As these are
// quite
// complicated
// indirections
- // we check for
+ // which one
+ // does not
+ // usually get
+ // right at
+ // first
+ // attempt we
+ // check for
// the internal
// consistency.
Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
ExcInternalError());
Assert (!neighbor_child->has_children(), ExcInternalError());
+ // We need to
+ // reset the
+ // ``un_v_matrix''
+ // on each
+ // subface
+ // because on
+ // each subface
+ // the ``un''
+ // belong to
+ // different
+ // neighboring
+ // cells.
+ un_v_matrix.clear();
+
// As already
// mentioned
// above for
- // this case
- // (case 2) we
- // employ the
- // FESubfaceValues
+ // the current
+ // case (case
+ // 2) we employ
+ // the
+ // ``FESubfaceValues''
// of the
// current
- // cell, here
+ // cell (here
// reinited for
// the current
// cell, face
- // and subface,
+ // and subface)
// and we
// employ the
// FEFaceValues
u_v_matrix,
un_v_matrix);
- // get dof
+ // Then we get
+ // the dof
// indices of
// the
// neighbor_child
// cell
neighbor_child->get_dof_indices (dofs_neighbor);
+ // and
// distribute
- // cell matrix
+ // ``un_v_matrix''
// to the
// system_matrix
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int k=0; k<dofs_per_cell; ++k)
system_matrix.add(dofs[i], dofs_neighbor[k],
un_v_matrix(i,k));
-
- // In the
- // ``assemble_face_term1''
- // function contributions to
- // the cell matrices and the
- // cell vector are only
- // ADDED. Therefore on each
- // subface we need to reset the
- // un_v_matrix
- // to zero, before assembling
- // the face terms corresponding
- // to the following neighbor_child cell.
- un_v_matrix.clear();
}
+ // End of ``if
+ // (face->has_children())''
}
- // End of ``if
- // (face->has_children())''
else
{
// We proceed with
// current cell.
if (neighbor->level() == cell->level())
{
- // Like before we
- // store which
- // number the
- // current cell has
- // in the list of
- // neighbors of the
- // neighboring
- // cell.
+ // Like before
+ // we store
+ // which number
+ // the current
+ // cell has in
+ // the list of
+ // neighbors of
+ // the
+ // neighboring
+ // cell.
const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
// We reinit
// the
- // FEFaceValues
+ // ``FEFaceValues''
// of the
// current and
// neighboring
// Reinit the
// appropriate
- // FEFaceValues
+ // ``FEFaceValues''
// and assemble
// the face
// terms.
un_v_matrix);
}
- // Get dof indices
- // of the
- // neighbor_child
+ // Now we get the
+ // dof indices of
+ // the
+ // ``neighbor_child''
// cell,
neighbor->get_dof_indices (dofs_neighbor);
- // distribute the
- // un_v_matrix,
+ // and distribute the
+ // ``un_v_matrix''.
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int k=0; k<dofs_per_cell; ++k)
system_matrix.add(dofs[i], dofs_neighbor[k],
un_v_matrix(i,k));
-
- // and clear the
- // ``un_v_matrix''
- // on each face.
- un_v_matrix.clear();
}
// End of ``face not at boundary'':
}
}
// Finally we distribute the
- // u_v_matrix,
+ // ``u_v_matrix''
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
- // the cell vector
+ // and the cell vector.
for (unsigned int i=0; i<dofs_per_cell; ++i)
right_hand_side(dofs[i]) += cell_vector(i);
-
- // and clear them both.
- u_v_matrix.clear ();
- cell_vector.clear ();
}
};
void DGMethod<dim>::assemble_system2 ()
{
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
- vector<unsigned int> dofs (dofs_per_cell);
- vector<unsigned int> dofs_neighbor (dofs_per_cell);
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
- UpdateFlags update_flags = UpdateFlags(update_values
- | update_gradients
- | update_q_points
- | update_JxW_values);
+ UpdateFlags update_flags = update_values
+ | update_gradients
+ | update_q_points
+ | update_JxW_values;
- UpdateFlags face_update_flags = UpdateFlags(update_values
- | update_q_points
- | update_JxW_values
- | update_normal_vectors);
-
- UpdateFlags neighbor_face_update_flags = UpdateFlags(update_values);
+ UpdateFlags face_update_flags = update_values
+ | update_q_points
+ | update_JxW_values
+ | update_normal_vectors;
+
+ UpdateFlags neighbor_face_update_flags = update_values;
// Here we do not need
// ``fe_v_face_neighbor'' as case 4
FullMatrix<double> u_v_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> un_v_matrix (dofs_per_cell, dofs_per_cell);
- // Additionally we need following
- // two cell matrices, both for face
- // term that include test function
- // ``vn'' (shape functions of the
- // neighboring cell). To be more
- // precise, the first matrix will
- // include the `u and vn terms' and
- // the second that will include the
- // `un and vn terms'.
+ // Additionally we need the
+ // following two cell matrices,
+ // both for face term that include
+ // test function ``vn'' (shape
+ // functions of the neighboring
+ // cell). To be more precise, the
+ // first matrix will include the `u
+ // and vn terms' and the second
+ // that will include the `un and vn
+ // terms'.
FullMatrix<double> u_vn_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> un_vn_matrix (dofs_per_cell, dofs_per_cell);
// The following lines are roughly
// the same as in the previous
// function.
- DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- DoFHandler<dim>::face_iterator face;
- DoFHandler<dim>::cell_iterator neighbor;
- DoFHandler<dim>::cell_iterator neighbor_child;
-
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
for (;cell!=endc; ++cell)
{
+ u_v_matrix.clear ();
+ cell_vector.clear ();
+
fe_v.reinit (cell);
dg.assemble_cell_term(fe_v,
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
- face = cell->face(face_no);
+ typename DoFHandler<dim>::face_iterator face=
+ cell->face(face_no);
// Case 1:
if (face->at_boundary())
else
{
Assert (cell->neighbor(face_no).state() == valid, ExcInternalError());
- neighbor = cell->neighbor(face_no);
-
+ typename DoFHandler<dim>::cell_iterator neighbor=
+ cell->neighbor(face_no);
// Case 2:
if (face->has_children())
{
subface_no<GeometryInfo<dim>::subfaces_per_face;
++subface_no)
{
- neighbor_child = neighbor->child(
- GeometryInfo<dim>::child_cell_on_face(neighbor2,subface_no));
+ typename DoFHandler<dim>::cell_iterator neighbor_child=
+ neighbor->child(GeometryInfo<dim>::child_cell_on_face(
+ neighbor2,subface_no));
Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
ExcInternalError());
Assert (!neighbor_child->has_children(), ExcInternalError());
-
+
+ un_v_matrix.clear();
+ u_vn_matrix.clear();
+ un_vn_matrix.clear();
+
fe_v_subface.reinit (cell, face_no, subface_no);
fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
un_vn_matrix(i,j));
}
-
- un_v_matrix.clear();
- u_vn_matrix.clear();
- un_vn_matrix.clear();
}
}
else
neighbor->index() > cell->index())
{
const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-
+
+ un_v_matrix.clear();
+ u_vn_matrix.clear();
+ un_vn_matrix.clear();
+
fe_v_face.reinit (cell, face_no);
fe_v_face_neighbor.reinit (neighbor, neighbor2);
system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
un_vn_matrix(i,j));
}
-
- un_v_matrix.clear();
- u_vn_matrix.clear();
- un_vn_matrix.clear();
}
// Due to rule b)
for (unsigned int i=0; i<dofs_per_cell; ++i)
right_hand_side(dofs[i]) += cell_vector(i);
-
- u_v_matrix.clear ();
- cell_vector.clear ();
}
};
- // For this simple solver we use the
+
+ // @sect3{All the rest}
+ //
+ // For this simple problem we use the
// simplest possible solver, called
- // richardson iteration, that
+ // Richardson iteration, that
// represents a simple defect
// correction. This, in combination
// with a block SSOR preconditioner,
// structur of system matrices
// arising from DG
// discretizations. The size of these
- // blocks are the number of DoFs
- // per cell. Here, we use a SSOR
+ // 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
// preconditioner step, it is wise
// to invert the diagonal blocks of
// the matrix before starting the
- // solver. Otherwise, the diagonal
- // blocks are inverted in each
+ // solver. Otherwise, it takes less
+ // memory, but the diagonal blocks
+ // are inverted in each
// preconditioner step,
// significantly slowing down the
// linear solving process.
// We refine the grid according to a
// very simple refinement criterion,
- // namely the gradients of the
- // solution. As here we consider the
- // DG(1) method (i.e. we use
- // piecewise bilinear shape
- // functions) we could simply compute
- // the gradients on each cell. But we
- // do not want to base our refinement
- // indicator on the gradients on each
- // cell only, but want to base them
- // also on jumps of the discontinuous
- // solution function over faces
- // between neighboring cells. The
- // simpliest way of doing that is to
- // compute approximative gradients by
+ // namely an approximation to the
+ // gradient of the solution. As here
+ // we consider the DG(1) method
+ // (i.e. we use piecewise bilinear
+ // shape functions) we could simply
+ // compute the gradients on each
+ // cell. But we do not want to base
+ // our refinement indicator on the
+ // gradients on each cell only, but
+ // want to base them also on jumps of
+ // the discontinuous solution
+ // function over faces between
+ // neighboring cells. The simpliest
+ // way of doing that is to compute
+ // approximative gradients by
// difference quotients including the
// cell under consideration and its
// neighbors. This is done by the
- // DerivativeApproximation class that
- // computes the approximate
+ // ``DerivativeApproximation'' class
+ // that computes the approximate
// gradients in a way similar to the
- // GradientEstimation described in
- // Step 9 of this tutorial. According
- // to the argumentation in Step 9,
- // here we consider
- // $h^{1+d/2}|\nabla_h
+ // ``GradientEstimation'' described
+ // in Step 9 of this tutorial. In
+ // fact, the
+ // ``DerivativeApproximation'' class
+ // was developed following the
+ // ``GradientEstimation'' class of
+ // Step 9. Relating to the
+ // discussion in Step 9, here we
+ // consider $h^{1+d/2}|\nabla_h
// u_h|$. Futhermore we note that we
- // do not consider approximate
- // second derivatives because
- // solutions to the linear advection
- // equation are in general not in H^2
- // but in H^1 (to be more precise, in
- // H^1_\beta) only.
+ // do not consider approximate second
+ // derivatives because solutions to
+ // the linear advection equation are
+ // in general not in H^2 but in H^1
+ // (to be more precise, in H^1_\beta)
+ // only.
template <int dim>
void DGMethod<dim>::refine_grid ()
{
- // The DerivativeApproximation
+ // The ``DerivativeApproximation''
// class computes the gradients to
// float precision. This is
// sufficient as they are
Assert (cycle < 10, ExcInternalError());
filename += ".gnuplot";
- cout << "Writing solution to <" << filename << ">..." << endl;
+ cout << "Writing solution to <" << filename << ">..." << endl << endl;
std::ofstream gnuplot_output (filename.c_str());
DataOut<dim> data_out;
// two solutions for equality.
solution1-=solution2;
const double difference=solution1.linfty_norm();
- if (difference<1e-13)
- cout << "solution1 and solution2 do not differ." << endl;
+ if (difference>1e-13)
+ cout << "solution1 and solution2 differ!!" << endl;
// Finally we perform the
// output.
}
}
-
-
+ // The following ``main'' function is
+ // similar to previous examples and
+ // need not to be commented on.
int main ()
{
- DGMethod<2> dgmethod_2d;
- dgmethod_2d.run ();
-
+ try
+ {
+ DGMethod<2> dgmethod;
+ dgmethod.run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+
return 0;
};
projects / dealii-svn.git / commitdiff