/* $Id$ */
-/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */
+/* Author: Wolfgang Bangerth, 1999, Guido Kanschat, 2011 */
/* $Id$ */
/* */
-/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010 by the deal.II authors */
+/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010, 2011 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
// suggest. Since they do not need to be
// called from outside, they are made
// private to this class.
+
private:
- void make_grid_and_dofs ();
+ void make_grid ();
+ void setup_system ();
void assemble_system ();
void solve ();
void output_results () const;
{}
- // @sect4{LaplaceProblem::make_grid_and_dofs}
+ // @sect4{LaplaceProblem::make_grid}
// Now, the first thing we've got to
// do is to generate the
// like to do our computation and
// number each vertex with a degree
// of freedom. We have seen this in
- // the previous examples before. Then
- // we have to set up space for the
- // system matrix and right hand side
- // of the discretized problem. This
- // is what this function does:
-void LaplaceProblem::make_grid_and_dofs ()
+ // the previous examples before.
+void LaplaceProblem::make_grid ()
{
// First create the grid and refine
// all cells five times. Since the
<< std::endl;
// Note the distinction between
// n_active_cells() and n_cells().
+}
+
+ // @sect4{LaplaceProblem::setup_system}
// Next we enumerate all the degrees of
- // freedom. This is done by using the
- // distribute_dofs function, as we have
+ // freedom and set up matrix and
+ // vector objects to hold the
+ // system data. Enumerating is done by using
+ // DoFHandler::distribute_dofs(), as we have
// seen in the step-2 example. Since we use
- // the <code>FE_Q</code> class with a polynomial
+ // the FE_Q class with a polynomial
// degree of 1, i.e. bilinear elements,
// this associates one degree of freedom
// with each vertex. While we're at
// generating output, let us also take a
// look at how many degrees of freedom are
// generated:
+void LaplaceProblem::setup_system ()
+{
dof_handler.distribute_dofs (fe);
std::cout << "Number of degrees of freedom: "
<< dof_handler.n_dofs()
// objects. That's too much, so there is one
// type of class that orchestrates
// information exchange between these three:
- // the <code>FEValues</code> class. If given one
+ // the FEValues class. If given one
// instance of each three of these objects,
// it will be able to provide you with
// information about values and gradients of
// actually need is given as a bitwise
// connection of flags as the third
// argument to the constructor of
- // <code>FEValues</code>. Since these values have to
+ // FEValues. Since these values have to
// be recomputed, or updated, every time we
// go to a new cell, all of these flags
// start with the prefix <code>update_</code> and
// then indicate what it actually is that
// we want updated. The flag to give if we
// want the values of the shape functions
- // computed is <code>update_values</code>; for the
+ // computed is #update_values; for the
// gradients it is
- // <code>update_gradients</code>. The determinants
+ // #update_gradients. The determinants
// of the Jacobians and the quadrature
// weights are always used together, so
// only the products (Jacobians times
// weights, or short <code>JxW</code>) are computed;
// since we need them, we have to list
- // <code>update_JxW_values</code> as well:
+ // #update_JxW_values as well:
FEValues<2> fe_values (fe, quadrature_formula,
update_values | update_gradients | update_JxW_values);
// The advantage of this proceeding is that
// determinant and the quadrature point
// weight (that one gets together by
// the call to
- // <code>fe_values.JxW</code>). Finally, this is
+ // FEValues::JxW() ). Finally, this is
// repeated for all shape functions
// phi_i and phi_j:
for (unsigned int i=0; i<dofs_per_cell; ++i)
// rather than projecting it onto the
// boundary. There is a function in the
// library which does exactly this:
- // <code>VectorTools::interpolate_boundary_values</code>. Its
+ // VectorTools::interpolate_boundary_values(). Its
// parameters are (omitting parameters for
// which default values exist and that we
// don't care about): the DoFHandler object
// the boundary.
//
// The function describing the boundary
- // values is an object of type <code>Function</code>
+ // values is an object of type Function
// or of a derived class. One of the
- // derived classes is <code>ZeroFunction</code>,
+ // derived classes is ZeroFunction,
// which describes (not unexpectedly) a
// function which is zero everywhere. We
// create such an object in-place and pass
- // it to the interpolate_boundary_values
+ // it to the VectorTools::interpolate_boundary_values()
// function.
//
// Finally, the output object is a
// First, we need to have an object that
// knows how to tell the CG algorithm when
// to stop. This is done by using a
- // <code>SolverControl</code> object, and as
+ // SolverControl object, and as
// stopping criterion we say: stop after a
// maximum of 1000 iterations (which is far
// more than is needed for 1089 variables;
// the one which stops the iteration:
SolverControl solver_control (1000, 1e-12);
// Then we need the solver itself. The
- // template parameters to the <code>SolverCG</code>
+ // template parameters to the SolverCG
// class are the matrix type and the type
// of the vectors, but the empty angle
// brackets indicate that we simply take
// To write the output to a file,
// we need an object which knows
// about output formats and the
- // like. This is the <code>DataOut</code> class,
+ // like. This is the DataOut class,
// and we need an object of that
// type:
DataOut<2> data_out;
// Now we have to tell it where to take the
// values from which it shall write. We
- // tell it which <code>DoFHandler</code> object to
+ // tell it which DoFHandler object to
// use, and the solution vector (and
// the name by which the solution variable
// shall appear in the output file). If
// handle. The reason is that we
// have separated the frontend
// (which knows about how to treat
- // <code>DoFHandler</code> objects and data
+ // DoFHandler objects and data
// vectors) from the back end (which
// knows many different output formats)
// and use an intermediate data
// to comment about:
void LaplaceProblem::run ()
{
- make_grid_and_dofs ();
+ make_grid ();
+ setup_system();
assemble_system ();
solve ();
output_results ();
// into the global namespace:
using namespace dealii;
- // @sect3{The <code>LaplaceProblem</code> class template}
+ // @sect3{The <code>Step4</code> class template}
// This is again the same
- // <code>LaplaceProblem</code> class as in the
+ // <code>Step4</code> class as in the
// previous example. The only
// difference is that we have now
// declared it as a class with a
// respectively. Apart from this,
// everything is as before.
template <int dim>
-class LaplaceProblem
+class Step4
{
public:
- LaplaceProblem ();
+ Step4 ();
void run ();
private:
- void make_grid_and_dofs ();
+ void make_grid ();
+ void setup_system();
void assemble_system ();
void solve ();
void output_results () const;
- // @sect3{Implementation of the <code>LaplaceProblem</code> class}
+ // @sect3{Implementation of the <code>Step4</code> class}
// Next for the implementation of the class
// template that makes use of the functions
// the time we define the template
// functions. Only later, the compiler will
// find a declaration of
- // <code>LaplaceProblem@<2@></code> (in the
+ // <code>Step4@<2@></code> (in the
// <code>main</code> function, actually) and
// compile the entire class with
// <code>dim</code> replaced by 2, a process
//
// In fact, the compiler will also find a
// declaration
- // <code>LaplaceProblem@<3@></code> in
+ // <code>Step4@<3@></code> in
// <code>main()</code>. This will cause it to
// again go back to the general
- // <code>LaplaceProblem@<dim@></code>
+ // <code>Step4@<dim@></code>
// template, replace all occurrences of
// <code>dim</code>, this time by 3, and
// compile the class a second time. Note that
// the two instantiations
- // <code>LaplaceProblem@<2@></code> and
- // <code>LaplaceProblem@<3@></code> are
+ // <code>Step4@<2@></code> and
+ // <code>Step4@<3@></code> are
// completely independent classes; their only
// common feature is that they are both
// instantiated from the same general
// completely independently).
- // @sect4{LaplaceProblem::LaplaceProblem}
+ // @sect4{Step4::Step4}
// After this introduction, here is the
- // constructor of the <code>LaplaceProblem</code>
+ // constructor of the <code>Step4</code>
// class. It specifies the desired polynomial
// degree of the finite elements and
// associates the DoFHandler to the
// triangulation just as in the previous
// example program, step-3:
template <int dim>
-LaplaceProblem<dim>::LaplaceProblem ()
+Step4<dim>::Step4 ()
:
fe (1),
dof_handler (triangulation)
{}
- // @sect4{LaplaceProblem::make_grid_and_dofs}
+ // @sect4{Step4::make_grid}
// Grid creation is something inherently
// dimension dependent. However, as long as
// solve on the square $[-1,1]\times [-1,1]$
// in 2D, or on the cube $[-1,1] \times
// [-1,1] \times [-1,1]$ in 3D; both can be
- // termed <code>hyper_cube</code>, so we may
+ // termed GridGenerator::hyper_cube(), so we may
// use the same function in whatever
// dimension we are. Of course, the functions
// that create a hypercube in two and three
// that is something you need not care
// about. Let the library handle the
// difficult things.
- //
- // Likewise, associating a degree of freedom
- // with each vertex is something which
- // certainly looks different in 2D and 3D,
- // but that does not need to bother you
- // either. This function therefore looks
- // exactly like in the previous example,
- // although it performs actions that in their
- // details are quite different if
- // <code>dim</code> happens to be 3. The only
- // significant difference from a user's
- // perspective is the number of cells
- // resulting, which is much higher in three
- // than in two space dimensions!
template <int dim>
-void LaplaceProblem<dim>::make_grid_and_dofs ()
+void Step4<dim>::make_grid ()
{
GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (4);
<< " Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
+}
+
+ // @sect4{Step4::setup_system}
+ // This function looks
+ // exactly like in the previous example,
+ // although it performs actions that in their
+ // details are quite different if
+ // <code>dim</code> happens to be 3. The only
+ // significant difference from a user's
+ // perspective is the number of cells
+ // resulting, which is much higher in three
+ // than in two space dimensions!
+template <int dim>
+void Step4<dim>::setup_system ()
+{
dof_handler.distribute_dofs (fe);
std::cout << " Number of degrees of freedom: "
}
- // @sect4{LaplaceProblem::assemble_system}
+ // @sect4{Step4::assemble_system}
// Unlike in the previous example, we
// would now like to use a
// don't have to care about most
// things.
template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
+void Step4<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(2);
// presently on (previously, we only
// required values and gradients of the
// shape function from the
- // <code>FEValues</code> object, as well as
+ // FEValues object, as well as
// the quadrature weights,
- // <code>JxW</code>). We can tell the
- // <code>FEValues</code> object to do for
+ // FEValues::JxW() ). We can tell the
+ // FEValues object to do for
// us by also giving it the
- // <code>update_quadrature_points</code>
+ // #update_quadrature_points
// flag:
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
// values in this example, unlike the one
// before. This is a simple task, we only
// have to replace the
- // <code>ZeroFunction</code> used there by
+ // ZeroFunction used there by
// an object of the class which describes
// the boundary values we would like to use
// (i.e. the <code>BoundaryValues</code>
}
- // @sect4{LaplaceProblem::solve}
+ // @sect4{Step4::solve}
// Solving the linear system of
// equations is something that looks
// function is copied verbatim from the
// previous example.
template <int dim>
-void LaplaceProblem<dim>::solve ()
+void Step4<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> solver (solver_control);
}
- // @sect4{LaplaceProblem::output_results}
+ // @sect4{Step4::output_results}
// This function also does what the
// respective one did in step-3. No changes
// than 2 or 3, but we neglect this here for
// the sake of brevity).
template <int dim>
-void LaplaceProblem<dim>::output_results () const
+void Step4<dim>::output_results () const
{
DataOut<dim> data_out;
- // @sect4{LaplaceProblem::run}
+ // @sect4{Step4::run}
// This is the function which has the
// top-level control over
// additional output, it is the same
// as for the previous example.
template <int dim>
-void LaplaceProblem<dim>::run ()
+void Step4<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
- make_grid_and_dofs();
+ make_grid();
+ setup_system ();
assemble_system ();
solve ();
output_results ();
// looks mostly like in step-3, but if you
// look at the code below, note how we first
// create a variable of type
- // <code>LaplaceProblem@<2@></code> (forcing
+ // <code>Step4@<2@></code> (forcing
// the compiler to compile the class template
// with <code>dim</code> replaced by
// <code>2</code>) and run a 2d simulation,
{
deallog.depth_console (0);
{
- LaplaceProblem<2> laplace_problem_2d;
+ Step4<2> laplace_problem_2d;
laplace_problem_2d.run ();
}
{
- LaplaceProblem<3> laplace_problem_3d;
+ Step4<3> laplace_problem_3d;
laplace_problem_3d.run ();
}