]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Do like in step-7: put everything into a namespace.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 19 Aug 2011 04:39:39 +0000 (04:39 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 19 Aug 2011 04:39:39 +0000 (04:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@24113 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-8/step-8.cc

index aecef3b5292dc5ea17b8cc7fce5b8f212dab3289..d617f803dfa649eea24ef99c43029ea1114dbd03 100644 (file)
@@ -3,7 +3,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010 by the deal.II authors */
+/*    Copyright (C) 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     */
 #include <iostream>
 
                                 // The last step is as in previous
-                                // programs:
-using namespace dealii;
-
-                                 // @sect3{The <code>ElasticProblem</code> class template}
-
-                                // The main class is, except for its
-                                // name, almost unchanged with
-                                // respect to the step-6 example.
-                                //
-                                // The only change is the use of a
-                                // different class for the <code>fe</code>
-                                // variable: Instead of a concrete
-                                // finite element class such as
-                                // <code>FE_Q</code>, we now use a more
-                                // generic one, <code>FESystem</code>. In
-                                // fact, <code>FESystem</code> is not really a
-                                // finite element itself in that it
-                                // does not implement shape functions
-                                // of its own.  Rather, it is a class
-                                // that can be used to stack several
-                                // other elements together to form
-                                // one vector-valued finite
-                                // element. In our case, we will
-                                // compose the vector-valued element
-                                // of <code>FE_Q(1)</code> objects, as shown
-                                // below in the constructor of this
-                                // class.
-template <int dim>
-class ElasticProblem
+                                // programs. In particular, just like in
+                                // step-7, we pack everything that's specific
+                                // to this program into a namespace of its
+                                // own.
+namespace Step8
 {
-  public:
-    ElasticProblem ();
-    ~ElasticProblem ();
-    void run ();
-
-  private:
-    void setup_system ();
-    void assemble_system ();
-    void solve ();
-    void refine_grid ();
-    void output_results (const unsigned int cycle) const;
-
-    Triangulation<dim>   triangulation;
-    DoFHandler<dim>      dof_handler;
-
-    FESystem<dim>        fe;
-
-    ConstraintMatrix     hanging_node_constraints;
-
-    SparsityPattern      sparsity_pattern;
-    SparseMatrix<double> system_matrix;
-
-    Vector<double>       solution;
-    Vector<double>       system_rhs;
-};
-
-
-                                 // @sect3{Right hand side values}
-
-                                // Before going over to the
-                                // implementation of the main class,
-                                // we declare and define the class
-                                // which describes the right hand
-                                // side. This time, the right hand
-                                // side is vector-valued, as is the
-                                // solution, so we will describe the
-                                // changes required for this in some
-                                // more detail.
-                                //
-                                // The first thing is that
-                                // vector-valued functions have to
-                                // have a constructor, since they
-                                // need to pass down to the base
-                                // class of how many components the
-                                // function consists. The default
-                                // value in the constructor of the
-                                // base class is one (i.e.: a scalar
-                                // function), which is why we did not
-                                // need not define a constructor for
-                                // the scalar function used in
-                                // previous programs.
-template <int dim>
-class RightHandSide :  public Function<dim>
-{
-  public:
-    RightHandSide ();
-
-                                    // The next change is that we
-                                    // want a replacement for the
-                                    // <code>value</code> function of the
-                                    // previous examples. There, a
-                                    // second parameter <code>component</code>
-                                    // was given, which denoted which
-                                    // component was requested. Here,
-                                    // we implement a function that
-                                    // returns the whole vector of
-                                    // values at the given place at
-                                    // once, in the second argument
-                                    // of the function. The obvious
-                                    // name for such a replacement
-                                    // function is <code>vector_value</code>.
-                                    //
-                                    // Secondly, in analogy to the
-                                    // <code>value_list</code> function, there
-                                    // is a function
-                                    // <code>vector_value_list</code>, which
-                                    // returns the values of the
-                                    // vector-valued function at
-                                    // several points at once:
-    virtual void vector_value (const Point<dim> &p,
-                              Vector<double>   &values) const;
-
-    virtual void vector_value_list (const std::vector<Point<dim> > &points,
-                                   std::vector<Vector<double> >   &value_list) const;
-};
-
-
-                                // This is the constructor of the
-                                // right hand side class. As said
-                                // above, it only passes down to the
-                                // base class the number of
-                                // components, which is <code>dim</code> in
-                                // the present case (one force
-                                // component in each of the <code>dim</code>
-                                // space directions).
-                                //
-                                // Some people would have moved the
-                                // definition of such a short
-                                // function right into the class
-                                // declaration. We do not do that, as
-                                // a matter of style: the deal.II
-                                // style guides require that class
-                                // declarations contain only
-                                // declarations, and that definitions
-                                // are always to be found
-                                // outside. This is, obviously, as
-                                // much as matter of taste as
-                                // indentation, but we try to be
-                                // consistent in this direction.
-template <int dim>
-RightHandSide<dim>::RightHandSide ()
-               :
-               Function<dim> (dim)
-{}
-
-
-                                // Next the function that returns
-                                // the whole vector of values at the
-                                // point <code>p</code> at once.
-                                //
-                                // To prevent cases where the return
-                                // vector has not previously been set
-                                // to the right size we test for this
-                                // case and otherwise throw an
-                                // exception at the beginning of the
-                                // function. Note that enforcing that
-                                // output arguments already have the
-                                // correct size is a convention in
-                                // deal.II, and enforced almost
-                                // everywhere. The reason is that we
-                                // would otherwise have to check at
-                                // the beginning of the function and
-                                // possibly change the size of the
-                                // output vector. This is expensive,
-                                // and would almost always be
-                                // unnecessary (the first call to the
-                                // function would set the vector to
-                                // the right size, and subsequent
-                                // calls would only have to do
-                                // redundant checks). In addition,
-                                // checking and possibly resizing the
-                                // vector is an operation that can
-                                // not be removed if we can't rely on
-                                // the assumption that the vector
-                                // already has the correct size; this
-                                // is in contract to the <code>Assert</code>
-                                // call that is completely removed if
-                                // the program is compiled in
-                                // optimized mode.
-                                //
-                                // Likewise, if by some accident
-                                // someone tried to compile and run
-                                // the program in only one space
-                                // dimension (in which the elastic
-                                // equations do not make much sense
-                                // since they reduce to the ordinary
-                                // Laplace equation), we terminate
-                                // the program in the second
-                                // assertion. The program will work
-                                // just fine in 3d, however.
-template <int dim>
-inline
-void RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                      Vector<double>   &values) const
-{
-  Assert (values.size() == dim,
-         ExcDimensionMismatch (values.size(), dim));
-  Assert (dim >= 2, ExcNotImplemented());
-
-                                  // The rest of the function
-                                  // implements computing force
-                                  // values. We will use a constant
-                                  // (unit) force in x-direction
-                                  // located in two little circles
-                                  // (or spheres, in 3d) around
-                                  // points (0.5,0) and (-0.5,0), and
-                                  // y-force in an area around the
-                                  // origin; in 3d, the z-component
-                                  // of these centers is zero as
-                                  // well.
-                                  //
-                                  // For this, let us first define
-                                  // two objects that denote the
-                                  // centers of these areas. Note
-                                  // that upon construction of the
-                                  // <code>Point</code> objects, all
-                                  // components are set to zero.
-  Point<dim> point_1, point_2;
-  point_1(0) = 0.5;
-  point_2(0) = -0.5;
-
-                                  // If now the point <code>p</code> is in a
-                                  // circle (sphere) of radius 0.2
-                                  // around one of these points, then
-                                  // set the force in x-direction to
-                                  // one, otherwise to zero:
-  if (((p-point_1).square() < 0.2*0.2) ||
-      ((p-point_2).square() < 0.2*0.2))
-    values(0) = 1;
-  else
-    values(0) = 0;
-
-                                  // Likewise, if <code>p</code> is in the
-                                  // vicinity of the origin, then set
-                                  // the y-force to 1, otherwise to
-                                  // zero:
-  if (p.square() < 0.2*0.2)
-    values(1) = 1;
-  else
-    values(1) = 0;
-}
+  using namespace dealii;
 
+                                  // @sect3{The <code>ElasticProblem</code> class template}
 
-
-                                // Now, this is the function of the
-                                // right hand side class that returns
-                                // the values at several points at
-                                // once. The function starts out with
-                                // checking that the number of input
-                                // and output arguments is equal (the
-                                // sizes of the individual output
-                                // vectors will be checked in the
-                                // function that we call further down
-                                // below). Next, we define an
-                                // abbreviation for the number of
-                                // points which we shall work on, to
-                                // make some things simpler below.
-template <int dim>
-void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
-                                           std::vector<Vector<double> >   &value_list) const
-{
-  Assert (value_list.size() == points.size(),
-         ExcDimensionMismatch (value_list.size(), points.size()));
-
-  const unsigned int n_points = points.size();
-
-                                  // Finally we treat each of the
-                                  // points. In one of the previous
-                                  // examples, we have explained why
-                                  // the
-                                  // <code>value_list</code>/<code>vector_value_list</code>
-                                  // function had been introduced: to
-                                  // prevent us from calling virtual
-                                  // functions too frequently. On the
-                                  // other hand, we now need to
-                                  // implement the same function
-                                  // twice, which can lead to
-                                  // confusion if one function is
-                                  // changed but the other is
-                                  // not.
+                                  // The main class is, except for its
+                                  // name, almost unchanged with
+                                  // respect to the step-6 example.
                                   //
-                                  // We can prevent this situation by
-                                  // calling
-                                  // <code>RightHandSide::vector_value</code>
-                                  // on each point in the input
-                                  // list. Note that by giving the
-                                  // full name of the function,
-                                  // including the class name, we
-                                  // instruct the compiler to
-                                  // explicitly call this function,
-                                  // and not to use the virtual
-                                  // function call mechanism that
-                                  // would be used if we had just
-                                  // called <code>vector_value</code>. This is
-                                  // important, since the compiler
-                                  // generally can't make any
-                                  // assumptions which function is
-                                  // called when using virtual
-                                  // functions, and it therefore
-                                  // can't inline the called function
-                                  // into the site of the call. On
-                                  // the contrary, here we give the
-                                  // fully qualified name, which
-                                  // bypasses the virtual function
-                                  // call, and consequently the
-                                  // compiler knows exactly which
-                                  // function is called and will
-                                  // inline above function into the
-                                  // present location. (Note that we
-                                  // have declared the
-                                  // <code>vector_value</code> function above
-                                  // <code>inline</code>, though modern
-                                  // compilers are also able to
-                                  // inline functions even if they
-                                  // have not been declared as
-                                  // inline).
+                                  // The only change is the use of a
+                                  // different class for the <code>fe</code>
+                                  // variable: Instead of a concrete
+                                  // finite element class such as
+                                  // <code>FE_Q</code>, we now use a more
+                                  // generic one, <code>FESystem</code>. In
+                                  // fact, <code>FESystem</code> is not really a
+                                  // finite element itself in that it
+                                  // does not implement shape functions
+                                  // of its own.  Rather, it is a class
+                                  // that can be used to stack several
+                                  // other elements together to form
+                                  // one vector-valued finite
+                                  // element. In our case, we will
+                                  // compose the vector-valued element
+                                  // of <code>FE_Q(1)</code> objects, as shown
+                                  // below in the constructor of this
+                                  // class.
+  template <int dim>
+  class ElasticProblem
+  {
+    public:
+      ElasticProblem ();
+      ~ElasticProblem ();
+      void run ();
+
+    private:
+      void setup_system ();
+      void assemble_system ();
+      void solve ();
+      void refine_grid ();
+      void output_results (const unsigned int cycle) const;
+
+      Triangulation<dim>   triangulation;
+      DoFHandler<dim>      dof_handler;
+
+      FESystem<dim>        fe;
+
+      ConstraintMatrix     hanging_node_constraints;
+
+      SparsityPattern      sparsity_pattern;
+      SparseMatrix<double> system_matrix;
+
+      Vector<double>       solution;
+      Vector<double>       system_rhs;
+  };
+
+
+                                  // @sect3{Right hand side values}
+
+                                  // Before going over to the
+                                  // implementation of the main class,
+                                  // we declare and define the class
+                                  // which describes the right hand
+                                  // side. This time, the right hand
+                                  // side is vector-valued, as is the
+                                  // solution, so we will describe the
+                                  // changes required for this in some
+                                  // more detail.
                                   //
-                                  // It is worth noting why we go to
-                                  // such length explaining what we
-                                  // do. Using this construct, we
-                                  // manage to avoid any
-                                  // inconsistency: if we want to
-                                  // change the right hand side
-                                  // function, it would be difficult
-                                  // to always remember that we
-                                  // always have to change two
-                                  // functions in the same way. Using
-                                  // this forwarding mechanism, we
-                                  // only have to change a single
-                                  // place (the <code>vector_value</code>
-                                  // function), and the second place
-                                  // (the <code>vector_value_list</code>
-                                  // function) will always be
-                                  // consistent with it. At the same
-                                  // time, using virtual function
-                                  // call bypassing, the code is no
-                                  // less efficient than if we had
-                                  // written it twice in the first
-                                  // place:
-  for (unsigned int p=0; p<n_points; ++p)
-    RightHandSide<dim>::vector_value (points[p],
-                                     value_list[p]);
-}
-
-
-
-                                 // @sect3{The <code>ElasticProblem</code> class implementation}
-
-                                 // @sect4{ElasticProblem::ElasticProblem}
-
-                                // Following is the constructor of
-                                // the main class. As said before, we
-                                // would like to construct a
-                                // vector-valued finite element that
-                                // is composed of several scalar
-                                // finite elements (i.e., we want to
-                                // build the vector-valued element so
-                                // that each of its vector components
-                                // consists of the shape functions of
-                                // a scalar element). Of course, the
-                                // number of scalar finite elements we
-                                // would like to stack together
-                                // equals the number of components
-                                // the solution function has, which
-                                // is <code>dim</code> since we consider
-                                // displacement in each space
-                                // direction. The <code>FESystem</code> class
-                                // can handle this: we pass it the
-                                // finite element of which we would
-                                // like to compose the system of, and
-                                // how often it shall be repeated:
-
-template <int dim>
-ElasticProblem<dim>::ElasticProblem ()
-               :
-               dof_handler (triangulation),
-               fe (FE_Q<dim>(1), dim)
-{}
-                                // In fact, the <code>FESystem</code> class
-                                // has several more constructors
-                                // which can perform more complex
-                                // operations than just stacking
-                                // together several scalar finite
-                                // elements of the same type into
-                                // one; we will get to know these
-                                // possibilities in later examples.
-
-
-
-                                 // @sect4{ElasticProblem::~ElasticProblem}
-
-                                // The destructor, on the other hand,
-                                // is exactly as in step-6:
-template <int dim>
-ElasticProblem<dim>::~ElasticProblem ()
-{
-  dof_handler.clear ();
-}
-
-
-                                 // @sect4{ElasticProblem::setup_system}
-
-                                // Setting up the system of equations
-                                // is identitical to the function
-                                // used in the step-6 example. The
-                                // <code>DoFHandler</code> class and all other
-                                // classes used here are fully aware
-                                // that the finite element we want to
-                                // use is vector-valued, and take
-                                // care of the vector-valuedness of
-                                // the finite element themselves. (In
-                                // fact, they do not, but this does
-                                // not need to bother you: since they
-                                // only need to know how many degrees
-                                // of freedom there are per vertex,
-                                // line and cell, and they do not ask
-                                // what they represent, i.e. whether
-                                // the finite element under
-                                // consideration is vector-valued or
-                                // whether it is, for example, a
-                                // scalar Hermite element with
-                                // several degrees of freedom on each
-                                // vertex).
-template <int dim>
-void ElasticProblem<dim>::setup_system ()
-{
-  dof_handler.distribute_dofs (fe);
-  hanging_node_constraints.clear ();
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                          hanging_node_constraints);
-  hanging_node_constraints.close ();
-  sparsity_pattern.reinit (dof_handler.n_dofs(),
-                          dof_handler.n_dofs(),
-                          dof_handler.max_couplings_between_dofs());
-  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
-
-  hanging_node_constraints.condense (sparsity_pattern);
-
-  sparsity_pattern.compress();
-
-  system_matrix.reinit (sparsity_pattern);
-
-  solution.reinit (dof_handler.n_dofs());
-  system_rhs.reinit (dof_handler.n_dofs());
-}
-
-
-                                 // @sect4{ElasticProblem::assemble_system}
-
-                                // The big changes in this program
-                                // are in the creation of matrix and
-                                // right hand side, since they are
-                                // problem-dependent. We will go
-                                // through that process step-by-step,
-                                // since it is a bit more complicated
-                                // than in previous examples.
-                                //
-                                // The first parts of this function
-                                // are the same as before, however:
-                                // setting up a suitable quadrature
-                                // formula, initializing an
-                                // <code>FEValues</code> object for the
-                                // (vector-valued) finite element we
-                                // use as well as the quadrature
-                                // object, and declaring a number of
-                                // auxiliary arrays. In addition, we
-                                // declare the ever same two
-                                // abbreviations: <code>n_q_points</code> and
-                                // <code>dofs_per_cell</code>. The number of
-                                // degrees of freedom per cell we now
-                                // obviously ask from the composed
-                                // finite element rather than from
-                                // the underlying scalar Q1
-                                // element. Here, it is <code>dim</code> times
-                                // the number of degrees of freedom
-                                // per cell of the Q1 element, though
-                                // this is not explicit knowledge we
-                                // need to care about:
-template <int dim>
-void ElasticProblem<dim>::assemble_system ()
-{
-  QGauss<dim>  quadrature_formula(2);
-
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                          update_values   | update_gradients |
-                           update_quadrature_points | update_JxW_values);
-
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
-
-  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double>       cell_rhs (dofs_per_cell);
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-                                  // As was shown in previous
-                                  // examples as well, we need a
-                                  // place where to store the values
-                                  // of the coefficients at all the
-                                  // quadrature points on a cell. In
-                                  // the present situation, we have
-                                  // two coefficients, lambda and mu.
-  std::vector<double>     lambda_values (n_q_points);
-  std::vector<double>     mu_values (n_q_points);
-
-                                  // Well, we could as well have
-                                  // omitted the above two arrays
-                                  // since we will use constant
-                                  // coefficients for both lambda and
-                                  // mu, which can be declared like
-                                  // this. They both represent
-                                  // functions always returning the
-                                  // constant value 1.0. Although we
-                                  // could omit the respective
-                                  // factors in the assemblage of the
-                                  // matrix, we use them here for
-                                  // purpose of demonstration.
-  ConstantFunction<dim> lambda(1.), mu(1.);
-
-                                  // Then again, we need to have the
-                                  // same for the right hand
-                                  // side. This is exactly as before
-                                  // in previous examples. However,
-                                  // we now have a vector-valued
-                                  // right hand side, which is why
-                                  // the data type of the
-                                  // <code>rhs_values</code> array is
-                                  // changed. We initialize it by
-                                  // <code>n_q_points</code> elements, each of
-                                  // which is a <code>Vector@<double@></code>
-                                  // with <code>dim</code> elements.
-  RightHandSide<dim>      right_hand_side;
-  std::vector<Vector<double> > rhs_values (n_q_points,
-                                          Vector<double>(dim));
-
-
-                                  // Now we can begin with the loop
-                                  // over all cells:
-  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                                endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      cell_matrix = 0;
-      cell_rhs = 0;
-
-      fe_values.reinit (cell);
-
-                                      // Next we get the values of
-                                      // the coefficients at the
-                                      // quadrature points. Likewise
-                                      // for the right hand side:
-      lambda.value_list (fe_values.get_quadrature_points(), lambda_values);
-      mu.value_list     (fe_values.get_quadrature_points(), mu_values);
-
-      right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
-                                        rhs_values);
-
-                                      // Then assemble the entries of
-                                      // the local stiffness matrix
-                                      // and right hand side
-                                      // vector. This follows almost
-                                      // one-to-one the pattern
-                                      // described in the
-                                      // introduction of this
-                                      // example.  One of the few
-                                      // comments in place is that we
-                                      // can compute the number
-                                      // <code>comp(i)</code>, i.e. the index
-                                      // of the only nonzero vector
-                                      // component of shape function
-                                      // <code>i</code> using the
-                                      // <code>fe.system_to_component_index(i).first</code>
-                                      // function call below.
+                                  // The first thing is that
+                                  // vector-valued functions have to
+                                  // have a constructor, since they
+                                  // need to pass down to the base
+                                  // class of how many components the
+                                  // function consists. The default
+                                  // value in the constructor of the
+                                  // base class is one (i.e.: a scalar
+                                  // function), which is why we did not
+                                  // need not define a constructor for
+                                  // the scalar function used in
+                                  // previous programs.
+  template <int dim>
+  class RightHandSide :  public Function<dim>
+  {
+    public:
+      RightHandSide ();
+
+                                      // The next change is that we
+                                      // want a replacement for the
+                                      // <code>value</code> function of the
+                                      // previous examples. There, a
+                                      // second parameter <code>component</code>
+                                      // was given, which denoted which
+                                      // component was requested. Here,
+                                      // we implement a function that
+                                      // returns the whole vector of
+                                      // values at the given place at
+                                      // once, in the second argument
+                                      // of the function. The obvious
+                                      // name for such a replacement
+                                      // function is <code>vector_value</code>.
                                       //
-                                      // (By accessing the
-                                      // <code>first</code> variable of
-                                      // the return value of the
-                                      // <code>system_to_component_index</code>
-                                      // function, you might
-                                      // already have guessed
-                                      // that there is more in
-                                      // it. In fact, the
-                                      // function returns a
-                                      // <code>std::pair@<unsigned int,
-                                      // unsigned int@></code>, of
-                                      // which the first element
-                                      // is <code>comp(i)</code> and the
-                                      // second is the value
-                                      // <code>base(i)</code> also noted
-                                      // in the introduction, i.e.
-                                      // the index
-                                      // of this shape function
-                                      // within all the shape
-                                      // functions that are nonzero
-                                      // in this component,
-                                      // i.e. <code>base(i)</code> in the
-                                      // diction of the
-                                      // introduction. This is not a
-                                      // number that we are usually
-                                      // interested in, however.)
-                                      //
-                                      // With this knowledge, we can
-                                      // assemble the local matrix
-                                      // contributions:
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       {
-         const unsigned int
-           component_i = fe.system_to_component_index(i).first;
-
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           {
-             const unsigned int
-               component_j = fe.system_to_component_index(j).first;
-
-             for (unsigned int q_point=0; q_point<n_q_points;
-                  ++q_point)
-               {
-                 cell_matrix(i,j)
-                   +=
-                                                    // The first term
-                                                    // is (lambda d_i
-                                                    // u_i, d_j v_j)
-                                                    // + (mu d_i u_j,
-                                                    // d_j v_i).
-                                                    // Note that
-                                                    // <code>shape_grad(i,q_point)</code>
-                                                    // returns the
-                                                    // gradient of
-                                                    // the only
-                                                    // nonzero
-                                                    // component of
-                                                    // the i-th shape
-                                                    // function at
-                                                    // quadrature
-                                                    // point
-                                                    // q_point. The
-                                                    // component
-                                                    // <code>comp(i)</code> of
-                                                    // the gradient,
-                                                    // which is the
-                                                    // derivative of
-                                                    // this only
-                                                    // nonzero vector
-                                                    // component of
-                                                    // the i-th shape
-                                                    // function with
-                                                    // respect to the
-                                                    // comp(i)th
-                                                    // coordinate is
-                                                    // accessed by
-                                                    // the appended
-                                                    // brackets.
-                   (
-                     (fe_values.shape_grad(i,q_point)[component_i] *
-                      fe_values.shape_grad(j,q_point)[component_j] *
-                      lambda_values[q_point])
-                     +
-                     (fe_values.shape_grad(i,q_point)[component_j] *
-                      fe_values.shape_grad(j,q_point)[component_i] *
-                      mu_values[q_point])
-                     +
-                                                      // The second term is
-                                                      // (mu nabla u_i, nabla v_j).
-                                                      // We need not
-                                                      // access a
-                                                      // specific
+                                      // Secondly, in analogy to the
+                                      // <code>value_list</code> function, there
+                                      // is a function
+                                      // <code>vector_value_list</code>, which
+                                      // returns the values of the
+                                      // vector-valued function at
+                                      // several points at once:
+      virtual void vector_value (const Point<dim> &p,
+                                Vector<double>   &values) const;
+
+      virtual void vector_value_list (const std::vector<Point<dim> > &points,
+                                     std::vector<Vector<double> >   &value_list) const;
+  };
+
+
+                                  // This is the constructor of the
+                                  // right hand side class. As said
+                                  // above, it only passes down to the
+                                  // base class the number of
+                                  // components, which is <code>dim</code> in
+                                  // the present case (one force
+                                  // component in each of the <code>dim</code>
+                                  // space directions).
+                                  //
+                                  // Some people would have moved the
+                                  // definition of such a short
+                                  // function right into the class
+                                  // declaration. We do not do that, as
+                                  // a matter of style: the deal.II
+                                  // style guides require that class
+                                  // declarations contain only
+                                  // declarations, and that definitions
+                                  // are always to be found
+                                  // outside. This is, obviously, as
+                                  // much as matter of taste as
+                                  // indentation, but we try to be
+                                  // consistent in this direction.
+  template <int dim>
+  RightHandSide<dim>::RightHandSide ()
+                 :
+                 Function<dim> (dim)
+  {}
+
+
+                                  // Next the function that returns
+                                  // the whole vector of values at the
+                                  // point <code>p</code> at once.
+                                  //
+                                  // To prevent cases where the return
+                                  // vector has not previously been set
+                                  // to the right size we test for this
+                                  // case and otherwise throw an
+                                  // exception at the beginning of the
+                                  // function. Note that enforcing that
+                                  // output arguments already have the
+                                  // correct size is a convention in
+                                  // deal.II, and enforced almost
+                                  // everywhere. The reason is that we
+                                  // would otherwise have to check at
+                                  // the beginning of the function and
+                                  // possibly change the size of the
+                                  // output vector. This is expensive,
+                                  // and would almost always be
+                                  // unnecessary (the first call to the
+                                  // function would set the vector to
+                                  // the right size, and subsequent
+                                  // calls would only have to do
+                                  // redundant checks). In addition,
+                                  // checking and possibly resizing the
+                                  // vector is an operation that can
+                                  // not be removed if we can't rely on
+                                  // the assumption that the vector
+                                  // already has the correct size; this
+                                  // is in contract to the <code>Assert</code>
+                                  // call that is completely removed if
+                                  // the program is compiled in
+                                  // optimized mode.
+                                  //
+                                  // Likewise, if by some accident
+                                  // someone tried to compile and run
+                                  // the program in only one space
+                                  // dimension (in which the elastic
+                                  // equations do not make much sense
+                                  // since they reduce to the ordinary
+                                  // Laplace equation), we terminate
+                                  // the program in the second
+                                  // assertion. The program will work
+                                  // just fine in 3d, however.
+  template <int dim>
+  inline
+  void RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                        Vector<double>   &values) const
+  {
+    Assert (values.size() == dim,
+           ExcDimensionMismatch (values.size(), dim));
+    Assert (dim >= 2, ExcNotImplemented());
+
+                                    // The rest of the function
+                                    // implements computing force
+                                    // values. We will use a constant
+                                    // (unit) force in x-direction
+                                    // located in two little circles
+                                    // (or spheres, in 3d) around
+                                    // points (0.5,0) and (-0.5,0), and
+                                    // y-force in an area around the
+                                    // origin; in 3d, the z-component
+                                    // of these centers is zero as
+                                    // well.
+                                    //
+                                    // For this, let us first define
+                                    // two objects that denote the
+                                    // centers of these areas. Note
+                                    // that upon construction of the
+                                    // <code>Point</code> objects, all
+                                    // components are set to zero.
+    Point<dim> point_1, point_2;
+    point_1(0) = 0.5;
+    point_2(0) = -0.5;
+
+                                    // If now the point <code>p</code> is in a
+                                    // circle (sphere) of radius 0.2
+                                    // around one of these points, then
+                                    // set the force in x-direction to
+                                    // one, otherwise to zero:
+    if (((p-point_1).square() < 0.2*0.2) ||
+       ((p-point_2).square() < 0.2*0.2))
+      values(0) = 1;
+    else
+      values(0) = 0;
+
+                                    // Likewise, if <code>p</code> is in the
+                                    // vicinity of the origin, then set
+                                    // the y-force to 1, otherwise to
+                                    // zero:
+    if (p.square() < 0.2*0.2)
+      values(1) = 1;
+    else
+      values(1) = 0;
+  }
+
+
+
+                                  // Now, this is the function of the
+                                  // right hand side class that returns
+                                  // the values at several points at
+                                  // once. The function starts out with
+                                  // checking that the number of input
+                                  // and output arguments is equal (the
+                                  // sizes of the individual output
+                                  // vectors will be checked in the
+                                  // function that we call further down
+                                  // below). Next, we define an
+                                  // abbreviation for the number of
+                                  // points which we shall work on, to
+                                  // make some things simpler below.
+  template <int dim>
+  void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                             std::vector<Vector<double> >   &value_list) const
+  {
+    Assert (value_list.size() == points.size(),
+           ExcDimensionMismatch (value_list.size(), points.size()));
+
+    const unsigned int n_points = points.size();
+
+                                    // Finally we treat each of the
+                                    // points. In one of the previous
+                                    // examples, we have explained why
+                                    // the
+                                    // <code>value_list</code>/<code>vector_value_list</code>
+                                    // function had been introduced: to
+                                    // prevent us from calling virtual
+                                    // functions too frequently. On the
+                                    // other hand, we now need to
+                                    // implement the same function
+                                    // twice, which can lead to
+                                    // confusion if one function is
+                                    // changed but the other is
+                                    // not.
+                                    //
+                                    // We can prevent this situation by
+                                    // calling
+                                    // <code>RightHandSide::vector_value</code>
+                                    // on each point in the input
+                                    // list. Note that by giving the
+                                    // full name of the function,
+                                    // including the class name, we
+                                    // instruct the compiler to
+                                    // explicitly call this function,
+                                    // and not to use the virtual
+                                    // function call mechanism that
+                                    // would be used if we had just
+                                    // called <code>vector_value</code>. This is
+                                    // important, since the compiler
+                                    // generally can't make any
+                                    // assumptions which function is
+                                    // called when using virtual
+                                    // functions, and it therefore
+                                    // can't inline the called function
+                                    // into the site of the call. On
+                                    // the contrary, here we give the
+                                    // fully qualified name, which
+                                    // bypasses the virtual function
+                                    // call, and consequently the
+                                    // compiler knows exactly which
+                                    // function is called and will
+                                    // inline above function into the
+                                    // present location. (Note that we
+                                    // have declared the
+                                    // <code>vector_value</code> function above
+                                    // <code>inline</code>, though modern
+                                    // compilers are also able to
+                                    // inline functions even if they
+                                    // have not been declared as
+                                    // inline).
+                                    //
+                                    // It is worth noting why we go to
+                                    // such length explaining what we
+                                    // do. Using this construct, we
+                                    // manage to avoid any
+                                    // inconsistency: if we want to
+                                    // change the right hand side
+                                    // function, it would be difficult
+                                    // to always remember that we
+                                    // always have to change two
+                                    // functions in the same way. Using
+                                    // this forwarding mechanism, we
+                                    // only have to change a single
+                                    // place (the <code>vector_value</code>
+                                    // function), and the second place
+                                    // (the <code>vector_value_list</code>
+                                    // function) will always be
+                                    // consistent with it. At the same
+                                    // time, using virtual function
+                                    // call bypassing, the code is no
+                                    // less efficient than if we had
+                                    // written it twice in the first
+                                    // place:
+    for (unsigned int p=0; p<n_points; ++p)
+      RightHandSide<dim>::vector_value (points[p],
+                                       value_list[p]);
+  }
+
+
+
+                                  // @sect3{The <code>ElasticProblem</code> class implementation}
+
+                                  // @sect4{ElasticProblem::ElasticProblem}
+
+                                  // Following is the constructor of
+                                  // the main class. As said before, we
+                                  // would like to construct a
+                                  // vector-valued finite element that
+                                  // is composed of several scalar
+                                  // finite elements (i.e., we want to
+                                  // build the vector-valued element so
+                                  // that each of its vector components
+                                  // consists of the shape functions of
+                                  // a scalar element). Of course, the
+                                  // number of scalar finite elements we
+                                  // would like to stack together
+                                  // equals the number of components
+                                  // the solution function has, which
+                                  // is <code>dim</code> since we consider
+                                  // displacement in each space
+                                  // direction. The <code>FESystem</code> class
+                                  // can handle this: we pass it the
+                                  // finite element of which we would
+                                  // like to compose the system of, and
+                                  // how often it shall be repeated:
+
+  template <int dim>
+  ElasticProblem<dim>::ElasticProblem ()
+                 :
+                 dof_handler (triangulation),
+                 fe (FE_Q<dim>(1), dim)
+  {}
+                                  // In fact, the <code>FESystem</code> class
+                                  // has several more constructors
+                                  // which can perform more complex
+                                  // operations than just stacking
+                                  // together several scalar finite
+                                  // elements of the same type into
+                                  // one; we will get to know these
+                                  // possibilities in later examples.
+
+
+
+                                  // @sect4{ElasticProblem::~ElasticProblem}
+
+                                  // The destructor, on the other hand,
+                                  // is exactly as in step-6:
+  template <int dim>
+  ElasticProblem<dim>::~ElasticProblem ()
+  {
+    dof_handler.clear ();
+  }
+
+
+                                  // @sect4{ElasticProblem::setup_system}
+
+                                  // Setting up the system of equations
+                                  // is identitical to the function
+                                  // used in the step-6 example. The
+                                  // <code>DoFHandler</code> class and all other
+                                  // classes used here are fully aware
+                                  // that the finite element we want to
+                                  // use is vector-valued, and take
+                                  // care of the vector-valuedness of
+                                  // the finite element themselves. (In
+                                  // fact, they do not, but this does
+                                  // not need to bother you: since they
+                                  // only need to know how many degrees
+                                  // of freedom there are per vertex,
+                                  // line and cell, and they do not ask
+                                  // what they represent, i.e. whether
+                                  // the finite element under
+                                  // consideration is vector-valued or
+                                  // whether it is, for example, a
+                                  // scalar Hermite element with
+                                  // several degrees of freedom on each
+                                  // vertex).
+  template <int dim>
+  void ElasticProblem<dim>::setup_system ()
+  {
+    dof_handler.distribute_dofs (fe);
+    hanging_node_constraints.clear ();
+    DoFTools::make_hanging_node_constraints (dof_handler,
+                                            hanging_node_constraints);
+    hanging_node_constraints.close ();
+    sparsity_pattern.reinit (dof_handler.n_dofs(),
+                            dof_handler.n_dofs(),
+                            dof_handler.max_couplings_between_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+    hanging_node_constraints.condense (sparsity_pattern);
+
+    sparsity_pattern.compress();
+
+    system_matrix.reinit (sparsity_pattern);
+
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+  }
+
+
+                                  // @sect4{ElasticProblem::assemble_system}
+
+                                  // The big changes in this program
+                                  // are in the creation of matrix and
+                                  // right hand side, since they are
+                                  // problem-dependent. We will go
+                                  // through that process step-by-step,
+                                  // since it is a bit more complicated
+                                  // than in previous examples.
+                                  //
+                                  // The first parts of this function
+                                  // are the same as before, however:
+                                  // setting up a suitable quadrature
+                                  // formula, initializing an
+                                  // <code>FEValues</code> object for the
+                                  // (vector-valued) finite element we
+                                  // use as well as the quadrature
+                                  // object, and declaring a number of
+                                  // auxiliary arrays. In addition, we
+                                  // declare the ever same two
+                                  // abbreviations: <code>n_q_points</code> and
+                                  // <code>dofs_per_cell</code>. The number of
+                                  // degrees of freedom per cell we now
+                                  // obviously ask from the composed
+                                  // finite element rather than from
+                                  // the underlying scalar Q1
+                                  // element. Here, it is <code>dim</code> times
+                                  // the number of degrees of freedom
+                                  // per cell of the Q1 element, though
+                                  // this is not explicit knowledge we
+                                  // need to care about:
+  template <int dim>
+  void ElasticProblem<dim>::assemble_system ()
+  {
+    QGauss<dim>  quadrature_formula(2);
+
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                            update_values   | update_gradients |
+                            update_quadrature_points | update_JxW_values);
+
+    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       cell_rhs (dofs_per_cell);
+
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+                                    // As was shown in previous
+                                    // examples as well, we need a
+                                    // place where to store the values
+                                    // of the coefficients at all the
+                                    // quadrature points on a cell. In
+                                    // the present situation, we have
+                                    // two coefficients, lambda and mu.
+    std::vector<double>     lambda_values (n_q_points);
+    std::vector<double>     mu_values (n_q_points);
+
+                                    // Well, we could as well have
+                                    // omitted the above two arrays
+                                    // since we will use constant
+                                    // coefficients for both lambda and
+                                    // mu, which can be declared like
+                                    // this. They both represent
+                                    // functions always returning the
+                                    // constant value 1.0. Although we
+                                    // could omit the respective
+                                    // factors in the assemblage of the
+                                    // matrix, we use them here for
+                                    // purpose of demonstration.
+    ConstantFunction<dim> lambda(1.), mu(1.);
+
+                                    // Then again, we need to have the
+                                    // same for the right hand
+                                    // side. This is exactly as before
+                                    // in previous examples. However,
+                                    // we now have a vector-valued
+                                    // right hand side, which is why
+                                    // the data type of the
+                                    // <code>rhs_values</code> array is
+                                    // changed. We initialize it by
+                                    // <code>n_q_points</code> elements, each of
+                                    // which is a <code>Vector@<double@></code>
+                                    // with <code>dim</code> elements.
+    RightHandSide<dim>      right_hand_side;
+    std::vector<Vector<double> > rhs_values (n_q_points,
+                                            Vector<double>(dim));
+
+
+                                    // Now we can begin with the loop
+                                    // over all cells:
+    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                  endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+       cell_matrix = 0;
+       cell_rhs = 0;
+
+       fe_values.reinit (cell);
+
+                                        // Next we get the values of
+                                        // the coefficients at the
+                                        // quadrature points. Likewise
+                                        // for the right hand side:
+       lambda.value_list (fe_values.get_quadrature_points(), lambda_values);
+       mu.value_list     (fe_values.get_quadrature_points(), mu_values);
+
+       right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
+                                          rhs_values);
+
+                                        // Then assemble the entries of
+                                        // the local stiffness matrix
+                                        // and right hand side
+                                        // vector. This follows almost
+                                        // one-to-one the pattern
+                                        // described in the
+                                        // introduction of this
+                                        // example.  One of the few
+                                        // comments in place is that we
+                                        // can compute the number
+                                        // <code>comp(i)</code>, i.e. the index
+                                        // of the only nonzero vector
+                                        // component of shape function
+                                        // <code>i</code> using the
+                                        // <code>fe.system_to_component_index(i).first</code>
+                                        // function call below.
+                                        //
+                                        // (By accessing the
+                                        // <code>first</code> variable of
+                                        // the return value of the
+                                        // <code>system_to_component_index</code>
+                                        // function, you might
+                                        // already have guessed
+                                        // that there is more in
+                                        // it. In fact, the
+                                        // function returns a
+                                        // <code>std::pair@<unsigned int,
+                                        // unsigned int@></code>, of
+                                        // which the first element
+                                        // is <code>comp(i)</code> and the
+                                        // second is the value
+                                        // <code>base(i)</code> also noted
+                                        // in the introduction, i.e.
+                                        // the index
+                                        // of this shape function
+                                        // within all the shape
+                                        // functions that are nonzero
+                                        // in this component,
+                                        // i.e. <code>base(i)</code> in the
+                                        // diction of the
+                                        // introduction. This is not a
+                                        // number that we are usually
+                                        // interested in, however.)
+                                        //
+                                        // With this knowledge, we can
+                                        // assemble the local matrix
+                                        // contributions:
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         {
+           const unsigned int
+             component_i = fe.system_to_component_index(i).first;
+
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             {
+               const unsigned int
+                 component_j = fe.system_to_component_index(j).first;
+
+               for (unsigned int q_point=0; q_point<n_q_points;
+                    ++q_point)
+                 {
+                   cell_matrix(i,j)
+                     +=
+                                                      // The first term
+                                                      // is (lambda d_i
+                                                      // u_i, d_j v_j)
+                                                      // + (mu d_i u_j,
+                                                      // d_j v_i).
+                                                      // Note that
+                                                      // <code>shape_grad(i,q_point)</code>
+                                                      // returns the
+                                                      // gradient of
+                                                      // the only
+                                                      // nonzero
                                                       // component of
-                                                      // the
-                                                      // gradient,
-                                                      // since we
-                                                      // only have to
-                                                      // compute the
-                                                      // scalar
-                                                      // product of
-                                                      // the two
-                                                      // gradients,
-                                                      // of which an
-                                                      // overloaded
-                                                      // version of
-                                                      // the
-                                                      // operator*
-                                                      // takes care,
-                                                      // as in
-                                                      // previous
-                                                      // examples.
-                                                      //
-                                                      // Note that by
-                                                      // using the ?:
-                                                      // operator, we
-                                                      // only do this
-                                                      // if comp(i)
-                                                      // equals
-                                                      // comp(j),
-                                                      // otherwise a
-                                                      // zero is
-                                                      // added (which
-                                                      // will be
-                                                      // optimized
-                                                      // away by the
-                                                      // compiler).
-                     ((component_i == component_j) ?
-                      (fe_values.shape_grad(i,q_point) *
-                       fe_values.shape_grad(j,q_point) *
-                       mu_values[q_point])  :
-                      0)
-                   )
-                   *
-                   fe_values.JxW(q_point);
-               }
-           }
-       }
-
-                                      // Assembling the right hand
-                                      // side is also just as
-                                      // discussed in the
-                                      // introduction:
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       {
-         const unsigned int
-           component_i = fe.system_to_component_index(i).first;
-
-         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-           cell_rhs(i) += fe_values.shape_value(i,q_point) *
-                          rhs_values[q_point](component_i) *
-                          fe_values.JxW(q_point);
-       }
-
-                                      // The transfer from local
-                                      // degrees of freedom into the
-                                      // global matrix and right hand
-                                      // side vector does not depend
-                                      // on the equation under
-                                      // consideration, and is thus
-                                      // the same as in all previous
-                                      // examples. The same holds for
-                                      // the elimination of hanging
-                                      // nodes from the matrix and
-                                      // right hand side, once we are
-                                      // done with assembling the
-                                      // entire linear system:
-      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],
-                              cell_matrix(i,j));
-
-         system_rhs(local_dof_indices[i]) += cell_rhs(i);
-       }
-    }
-
-  hanging_node_constraints.condense (system_matrix);
-  hanging_node_constraints.condense (system_rhs);
-
-                                  // The interpolation of the
-                                  // boundary values needs a small
-                                  // modification: since the solution
-                                  // function is vector-valued, so
-                                  // need to be the boundary
-                                  // values. The <code>ZeroFunction</code>
-                                  // constructor accepts a parameter
-                                  // that tells it that it shall
-                                  // represent a vector valued,
-                                  // constant zero function with that
-                                  // many components. By default,
-                                  // this parameter is equal to one,
-                                  // in which case the
-                                  // <code>ZeroFunction</code> object would
-                                  // represent a scalar
-                                  // function. Since the solution
-                                  // vector has <code>dim</code> components,
-                                  // we need to pass <code>dim</code> as
-                                  // number of components to the zero
-                                  // function as well.
-  std::map<unsigned int,double> boundary_values;
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(dim),
-                                           boundary_values);
-  MatrixTools::apply_boundary_values (boundary_values,
-                                     system_matrix,
-                                     solution,
-                                     system_rhs);
-}
-
-
-
-                                 // @sect4{ElasticProblem::solve}
-
-                                // The solver does not care about
-                                // where the system of equations
-                                // comes, as long as it stays
-                                // positive definite and symmetric
-                                // (which are the requirements for
-                                // the use of the CG solver), which
-                                // the system indeed is. Therefore,
-                                // we need not change anything.
-template <int dim>
-void ElasticProblem<dim>::solve ()
-{
-  SolverControl           solver_control (1000, 1e-12);
-  SolverCG<>              cg (solver_control);
-
-  PreconditionSSOR<> preconditioner;
-  preconditioner.initialize(system_matrix, 1.2);
-
-  cg.solve (system_matrix, solution, system_rhs,
-           preconditioner);
-
-  hanging_node_constraints.distribute (solution);
-}
-
-
-                                 // @sect4{ElasticProblem::refine_grid}
-
-                                // The function that does the
-                                // refinement of the grid is the same
-                                // as in the step-6 example. The
-                                // quadrature formula is adapted to
-                                // the linear elements again. Note
-                                // that the error estimator by
-                                // default adds up the estimated
-                                // obtained from all components of
-                                // the finite element solution, i.e.,
-                                // it uses the displacement in all
-                                // directions with the same
-                                // weight. If we would like the grid
-                                // to be adapted to the
-                                // x-displacement only, we could pass
-                                // the function an additional
-                                // parameter which tells it to do so
-                                // and do not consider the
-                                // displacements in all other
-                                // directions for the error
-                                // indicators. However, for the
-                                // current problem, it seems
-                                // appropriate to consider all
-                                // displacement components with equal
-                                // weight.
-template <int dim>
-void ElasticProblem<dim>::refine_grid ()
-{
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-  typename FunctionMap<dim>::type neumann_boundary;
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                     QGauss<dim-1>(2),
-                                     neumann_boundary,
-                                     solution,
-                                     estimated_error_per_cell);
-
-  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                  estimated_error_per_cell,
-                                                  0.3, 0.03);
-
-  triangulation.execute_coarsening_and_refinement ();
-}
-
-
-                                 // @sect4{ElasticProblem::output_results}
-
-                                // The output happens mostly as has
-                                // been shown in previous examples
-                                // already. The only difference is
-                                // that the solution function is
-                                // vector valued. The <code>DataOut</code>
-                                // class takes care of this
-                                // automatically, but we have to give
-                                // each component of the solution
-                                // vector a different name.
-template <int dim>
-void ElasticProblem<dim>::output_results (const unsigned int cycle) const
-{
-  std::string filename = "solution-";
-  filename += ('0' + cycle);
-  Assert (cycle < 10, ExcInternalError());
-
-  filename += ".gmv";
-  std::ofstream output (filename.c_str());
-
-  DataOut<dim> data_out;
-  data_out.attach_dof_handler (dof_handler);
-
-
-
-                                  // As said above, we need a
-                                  // different name for each
-                                  // component of the solution
-                                  // function. To pass one name for
-                                  // each component, a vector of
-                                  // strings is used. Since the
-                                  // number of components is the same
-                                  // as the number of dimensions we
-                                  // are working in, the following
-                                  // <code>switch</code> statement is used.
+                                                      // the i-th shape
+                                                      // function at
+                                                      // quadrature
+                                                      // point
+                                                      // q_point. The
+                                                      // component
+                                                      // <code>comp(i)</code> of
+                                                      // the gradient,
+                                                      // which is the
+                                                      // derivative of
+                                                      // this only
+                                                      // nonzero vector
+                                                      // component of
+                                                      // the i-th shape
+                                                      // function with
+                                                      // respect to the
+                                                      // comp(i)th
+                                                      // coordinate is
+                                                      // accessed by
+                                                      // the appended
+                                                      // brackets.
+                     (
+                       (fe_values.shape_grad(i,q_point)[component_i] *
+                        fe_values.shape_grad(j,q_point)[component_j] *
+                        lambda_values[q_point])
+                       +
+                       (fe_values.shape_grad(i,q_point)[component_j] *
+                        fe_values.shape_grad(j,q_point)[component_i] *
+                        mu_values[q_point])
+                       +
+                                                        // The second term is
+                                                        // (mu nabla u_i, nabla v_j).
+                                                        // We need not
+                                                        // access a
+                                                        // specific
+                                                        // component of
+                                                        // the
+                                                        // gradient,
+                                                        // since we
+                                                        // only have to
+                                                        // compute the
+                                                        // scalar
+                                                        // product of
+                                                        // the two
+                                                        // gradients,
+                                                        // of which an
+                                                        // overloaded
+                                                        // version of
+                                                        // the
+                                                        // operator*
+                                                        // takes care,
+                                                        // as in
+                                                        // previous
+                                                        // examples.
+                                                        //
+                                                        // Note that by
+                                                        // using the ?:
+                                                        // operator, we
+                                                        // only do this
+                                                        // if comp(i)
+                                                        // equals
+                                                        // comp(j),
+                                                        // otherwise a
+                                                        // zero is
+                                                        // added (which
+                                                        // will be
+                                                        // optimized
+                                                        // away by the
+                                                        // compiler).
+                       ((component_i == component_j) ?
+                        (fe_values.shape_grad(i,q_point) *
+                         fe_values.shape_grad(j,q_point) *
+                         mu_values[q_point])  :
+                        0)
+                     )
+                     *
+                     fe_values.JxW(q_point);
+                 }
+             }
+         }
+
+                                        // Assembling the right hand
+                                        // side is also just as
+                                        // discussed in the
+                                        // introduction:
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         {
+           const unsigned int
+             component_i = fe.system_to_component_index(i).first;
+
+           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+             cell_rhs(i) += fe_values.shape_value(i,q_point) *
+                            rhs_values[q_point](component_i) *
+                            fe_values.JxW(q_point);
+         }
+
+                                        // The transfer from local
+                                        // degrees of freedom into the
+                                        // global matrix and right hand
+                                        // side vector does not depend
+                                        // on the equation under
+                                        // consideration, and is thus
+                                        // the same as in all previous
+                                        // examples. The same holds for
+                                        // the elimination of hanging
+                                        // nodes from the matrix and
+                                        // right hand side, once we are
+                                        // done with assembling the
+                                        // entire linear system:
+       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],
+                                cell_matrix(i,j));
+
+           system_rhs(local_dof_indices[i]) += cell_rhs(i);
+         }
+      }
+
+    hanging_node_constraints.condense (system_matrix);
+    hanging_node_constraints.condense (system_rhs);
+
+                                    // The interpolation of the
+                                    // boundary values needs a small
+                                    // modification: since the solution
+                                    // function is vector-valued, so
+                                    // need to be the boundary
+                                    // values. The <code>ZeroFunction</code>
+                                    // constructor accepts a parameter
+                                    // that tells it that it shall
+                                    // represent a vector valued,
+                                    // constant zero function with that
+                                    // many components. By default,
+                                    // this parameter is equal to one,
+                                    // in which case the
+                                    // <code>ZeroFunction</code> object would
+                                    // represent a scalar
+                                    // function. Since the solution
+                                    // vector has <code>dim</code> components,
+                                    // we need to pass <code>dim</code> as
+                                    // number of components to the zero
+                                    // function as well.
+    std::map<unsigned int,double> boundary_values;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             ZeroFunction<dim>(dim),
+                                             boundary_values);
+    MatrixTools::apply_boundary_values (boundary_values,
+                                       system_matrix,
+                                       solution,
+                                       system_rhs);
+  }
+
+
+
+                                  // @sect4{ElasticProblem::solve}
+
+                                  // The solver does not care about
+                                  // where the system of equations
+                                  // comes, as long as it stays
+                                  // positive definite and symmetric
+                                  // (which are the requirements for
+                                  // the use of the CG solver), which
+                                  // the system indeed is. Therefore,
+                                  // we need not change anything.
+  template <int dim>
+  void ElasticProblem<dim>::solve ()
+  {
+    SolverControl           solver_control (1000, 1e-12);
+    SolverCG<>              cg (solver_control);
+
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+
+    cg.solve (system_matrix, solution, system_rhs,
+             preconditioner);
+
+    hanging_node_constraints.distribute (solution);
+  }
+
+
+                                  // @sect4{ElasticProblem::refine_grid}
+
+                                  // The function that does the
+                                  // refinement of the grid is the same
+                                  // as in the step-6 example. The
+                                  // quadrature formula is adapted to
+                                  // the linear elements again. Note
+                                  // that the error estimator by
+                                  // default adds up the estimated
+                                  // obtained from all components of
+                                  // the finite element solution, i.e.,
+                                  // it uses the displacement in all
+                                  // directions with the same
+                                  // weight. If we would like the grid
+                                  // to be adapted to the
+                                  // x-displacement only, we could pass
+                                  // the function an additional
+                                  // parameter which tells it to do so
+                                  // and do not consider the
+                                  // displacements in all other
+                                  // directions for the error
+                                  // indicators. However, for the
+                                  // current problem, it seems
+                                  // appropriate to consider all
+                                  // displacement components with equal
+                                  // weight.
+  template <int dim>
+  void ElasticProblem<dim>::refine_grid ()
+  {
+    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+    typename FunctionMap<dim>::type neumann_boundary;
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                       QGauss<dim-1>(2),
+                                       neumann_boundary,
+                                       solution,
+                                       estimated_error_per_cell);
+
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.3, 0.03);
+
+    triangulation.execute_coarsening_and_refinement ();
+  }
+
+
+                                  // @sect4{ElasticProblem::output_results}
+
+                                  // The output happens mostly as has
+                                  // been shown in previous examples
+                                  // already. The only difference is
+                                  // that the solution function is
+                                  // vector valued. The <code>DataOut</code>
+                                  // class takes care of this
+                                  // automatically, but we have to give
+                                  // each component of the solution
+                                  // vector a different name.
+  template <int dim>
+  void ElasticProblem<dim>::output_results (const unsigned int cycle) const
+  {
+    std::string filename = "solution-";
+    filename += ('0' + cycle);
+    Assert (cycle < 10, ExcInternalError());
+
+    filename += ".gmv";
+    std::ofstream output (filename.c_str());
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+
+
+
+                                    // As said above, we need a
+                                    // different name for each
+                                    // component of the solution
+                                    // function. To pass one name for
+                                    // each component, a vector of
+                                    // strings is used. Since the
+                                    // number of components is the same
+                                    // as the number of dimensions we
+                                    // are working in, the following
+                                    // <code>switch</code> statement is used.
+                                    //
+                                    // We note that some graphics
+                                    // programs have restriction as to
+                                    // what characters are allowed in
+                                    // the names of variables. The
+                                    // library therefore supports only
+                                    // the minimal subset of these
+                                    // characters that is supported by
+                                    // all programs. Basically, these
+                                    // are letters, numbers,
+                                    // underscores, and some other
+                                    // characters, but in particular no
+                                    // whitespace and minus/hyphen. The
+                                    // library will throw an exception
+                                    // otherwise, at least if in debug
+                                    // mode.
+                                    //
+                                    // After listing the 1d, 2d, and 3d
+                                    // case, it is good style to let
+                                    // the program die if we run upon a
+                                    // case which we did not
+                                    // consider. Remember that the
+                                    // <code>Assert</code> macro generates an
+                                    // exception if the condition in
+                                    // the first parameter is not
+                                    // satisfied. Of course, the
+                                    // condition <code>false</code> can never be
+                                    // satisfied, so the program will
+                                    // always abort whenever it gets to
+                                    // the default statement:
+    std::vector<std::string> solution_names;
+    switch (dim)
+      {
+       case 1:
+             solution_names.push_back ("displacement");
+             break;
+       case 2:
+             solution_names.push_back ("x_displacement");
+             solution_names.push_back ("y_displacement");
+             break;
+       case 3:
+             solution_names.push_back ("x_displacement");
+             solution_names.push_back ("y_displacement");
+             solution_names.push_back ("z_displacement");
+             break;
+       default:
+             Assert (false, ExcNotImplemented());
+      }
+
+                                    // After setting up the names for
+                                    // the different components of the
+                                    // solution vector, we can add the
+                                    // solution vector to the list of
+                                    // data vectors scheduled for
+                                    // output. Note that the following
+                                    // function takes a vector of
+                                    // strings as second argument,
+                                    // whereas the one which we have
+                                    // used in all previous examples
+                                    // accepted a string there. In
+                                    // fact, the latter function is
+                                    // only a shortcut for the function
+                                    // which we call here: it puts the
+                                    // single string that is passed to
+                                    // it into a vector of strings with
+                                    // only one element and forwards
+                                    // that to the other function.
+    data_out.add_data_vector (solution, solution_names);
+    data_out.build_patches ();
+    data_out.write_gmv (output);
+  }
+
+
+
+                                  // @sect4{ElasticProblem::run}
+
+                                  // The <code>run</code> function does the same
+                                  // things as in step-6, for
+                                  // example. This time, we use the
+                                  // square [-1,1]^d as domain, and we
+                                  // refine it twice globally before
+                                  // starting the first iteration.
                                   //
-                                  // We note that some graphics
-                                  // programs have restriction as to
-                                  // what characters are allowed in
-                                  // the names of variables. The
-                                  // library therefore supports only
-                                  // the minimal subset of these
-                                  // characters that is supported by
-                                  // all programs. Basically, these
-                                  // are letters, numbers,
-                                  // underscores, and some other
-                                  // characters, but in particular no
-                                  // whitespace and minus/hyphen. The
-                                  // library will throw an exception
-                                  // otherwise, at least if in debug
-                                  // mode.
+                                  // The reason is the following: we
+                                  // use the <code>Gauss</code> quadrature
+                                  // formula with two points in each
+                                  // direction for integration of the
+                                  // right hand side; that means that
+                                  // there are four quadrature points
+                                  // on each cell (in 2D). If we only
+                                  // refine the initial grid once
+                                  // globally, then there will be only
+                                  // four quadrature points in each
+                                  // direction on the domain. However,
+                                  // the right hand side function was
+                                  // chosen to be rather localized and
+                                  // in that case all quadrature points
+                                  // lie outside the support of the
+                                  // right hand side function. The
+                                  // right hand side vector will then
+                                  // contain only zeroes and the
+                                  // solution of the system of
+                                  // equations is the zero vector,
+                                  // i.e. a finite element function
+                                  // that it zero everywhere. We should
+                                  // not be surprised about such things
+                                  // happening, since we have chosen an
+                                  // initial grid that is totally
+                                  // unsuitable for the problem at
+                                  // hand.
                                   //
-                                  // After listing the 1d, 2d, and 3d
-                                  // case, it is good style to let
-                                  // the program die if we run upon a
-                                  // case which we did not
-                                  // consider. Remember that the
-                                  // <code>Assert</code> macro generates an
-                                  // exception if the condition in
-                                  // the first parameter is not
-                                  // satisfied. Of course, the
-                                  // condition <code>false</code> can never be
-                                  // satisfied, so the program will
-                                  // always abort whenever it gets to
-                                  // the default statement:
-  std::vector<std::string> solution_names;
-  switch (dim)
-    {
-      case 1:
-           solution_names.push_back ("displacement");
-           break;
-      case 2:
-           solution_names.push_back ("x_displacement");
-           solution_names.push_back ("y_displacement");
-           break;
-      case 3:
-           solution_names.push_back ("x_displacement");
-           solution_names.push_back ("y_displacement");
-           solution_names.push_back ("z_displacement");
-           break;
-      default:
-           Assert (false, ExcNotImplemented());
-    }
-
-                                  // After setting up the names for
-                                  // the different components of the
-                                  // solution vector, we can add the
-                                  // solution vector to the list of
-                                  // data vectors scheduled for
-                                  // output. Note that the following
-                                  // function takes a vector of
-                                  // strings as second argument,
-                                  // whereas the one which we have
-                                  // used in all previous examples
-                                  // accepted a string there. In
-                                  // fact, the latter function is
-                                  // only a shortcut for the function
-                                  // which we call here: it puts the
-                                  // single string that is passed to
-                                  // it into a vector of strings with
-                                  // only one element and forwards
-                                  // that to the other function.
-  data_out.add_data_vector (solution, solution_names);
-  data_out.build_patches ();
-  data_out.write_gmv (output);
-}
-
-
-
-                                 // @sect4{ElasticProblem::run}
-
-                                // The <code>run</code> function does the same
-                                // things as in step-6, for
-                                // example. This time, we use the
-                                // square [-1,1]^d as domain, and we
-                                // refine it twice globally before
-                                // starting the first iteration.
-                                //
-                                // The reason is the following: we
-                                // use the <code>Gauss</code> quadrature
-                                // formula with two points in each
-                                // direction for integration of the
-                                // right hand side; that means that
-                                // there are four quadrature points
-                                // on each cell (in 2D). If we only
-                                // refine the initial grid once
-                                // globally, then there will be only
-                                // four quadrature points in each
-                                // direction on the domain. However,
-                                // the right hand side function was
-                                // chosen to be rather localized and
-                                // in that case all quadrature points
-                                // lie outside the support of the
-                                // right hand side function. The
-                                // right hand side vector will then
-                                // contain only zeroes and the
-                                // solution of the system of
-                                // equations is the zero vector,
-                                // i.e. a finite element function
-                                // that it zero everywhere. We should
-                                // not be surprised about such things
-                                // happening, since we have chosen an
-                                // initial grid that is totally
-                                // unsuitable for the problem at
-                                // hand.
-                                //
-                                // The unfortunate thing is that if
-                                // the discrete solution is constant,
-                                // then the error indicators computed
-                                // by the <code>KellyErrorEstimator</code>
-                                // class are zero for each cell as
-                                // well, and the call to
-                                // <code>refine_and_coarsen_fixed_number</code>
-                                // on the <code>triangulation</code> object
-                                // will not flag any cells for
-                                // refinement (why should it if the
-                                // indicated error is zero for each
-                                // cell?). The grid in the next
-                                // iteration will therefore consist
-                                // of four cells only as well, and
-                                // the same problem occurs again.
-                                //
-                                // The conclusion needs to be: while
-                                // of course we will not choose the
-                                // initial grid to be well-suited for
-                                // the accurate solution of the
-                                // problem, we must at least choose
-                                // it such that it has the chance to
-                                // capture the most striking features
-                                // of the solution. In this case, it
-                                // needs to be able to see the right
-                                // hand side. Thus, we refine twice
-                                // globally. (Note that the
-                                // <code>refine_global</code> function is not
-                                // part of the <code>GridRefinement</code>
-                                // class in which
-                                // <code>refine_and_coarsen_fixed_number</code>
-                                // is declared, for example. The
-                                // reason is first that it is not an
-                                // algorithm that computed refinement
-                                // flags from indicators, but more
-                                // importantly that it actually
-                                // performs the refinement, in
-                                // contrast to the functions in
-                                // <code>GridRefinement</code> that only flag
-                                // cells without actually refining
-                                // the grid.)
-template <int dim>
-void ElasticProblem<dim>::run ()
-{
-  for (unsigned int cycle=0; cycle<8; ++cycle)
-    {
-      std::cout << "Cycle " << cycle << ':' << std::endl;
-
-      if (cycle == 0)
-       {
-         GridGenerator::hyper_cube (triangulation, -1, 1);
-         triangulation.refine_global (2);
-       }
-      else
-       refine_grid ();
-
-      std::cout << "   Number of active cells:       "
-               << triangulation.n_active_cells()
-               << std::endl;
-
-      setup_system ();
-
-      std::cout << "   Number of degrees of freedom: "
-               << dof_handler.n_dofs()
-               << std::endl;
-
-      assemble_system ();
-      solve ();
-      output_results (cycle);
-    }
+                                  // The unfortunate thing is that if
+                                  // the discrete solution is constant,
+                                  // then the error indicators computed
+                                  // by the <code>KellyErrorEstimator</code>
+                                  // class are zero for each cell as
+                                  // well, and the call to
+                                  // <code>refine_and_coarsen_fixed_number</code>
+                                  // on the <code>triangulation</code> object
+                                  // will not flag any cells for
+                                  // refinement (why should it if the
+                                  // indicated error is zero for each
+                                  // cell?). The grid in the next
+                                  // iteration will therefore consist
+                                  // of four cells only as well, and
+                                  // the same problem occurs again.
+                                  //
+                                  // The conclusion needs to be: while
+                                  // of course we will not choose the
+                                  // initial grid to be well-suited for
+                                  // the accurate solution of the
+                                  // problem, we must at least choose
+                                  // it such that it has the chance to
+                                  // capture the most striking features
+                                  // of the solution. In this case, it
+                                  // needs to be able to see the right
+                                  // hand side. Thus, we refine twice
+                                  // globally. (Note that the
+                                  // <code>refine_global</code> function is not
+                                  // part of the <code>GridRefinement</code>
+                                  // class in which
+                                  // <code>refine_and_coarsen_fixed_number</code>
+                                  // is declared, for example. The
+                                  // reason is first that it is not an
+                                  // algorithm that computed refinement
+                                  // flags from indicators, but more
+                                  // importantly that it actually
+                                  // performs the refinement, in
+                                  // contrast to the functions in
+                                  // <code>GridRefinement</code> that only flag
+                                  // cells without actually refining
+                                  // the grid.)
+  template <int dim>
+  void ElasticProblem<dim>::run ()
+  {
+    for (unsigned int cycle=0; cycle<8; ++cycle)
+      {
+       std::cout << "Cycle " << cycle << ':' << std::endl;
+
+       if (cycle == 0)
+         {
+           GridGenerator::hyper_cube (triangulation, -1, 1);
+           triangulation.refine_global (2);
+         }
+       else
+         refine_grid ();
+
+       std::cout << "   Number of active cells:       "
+                 << triangulation.n_active_cells()
+                 << std::endl;
+
+       setup_system ();
+
+       std::cout << "   Number of degrees of freedom: "
+                 << dof_handler.n_dofs()
+                 << std::endl;
+
+       assemble_system ();
+       solve ();
+       output_results (cycle);
+      }
+  }
 }
 
                                  // @sect3{The <code>main</code> function}
 
-                                // The main function is again exactly
-                                // like in step-6 (apart from the
-                                // changed class names, of course).
+                                // After closing the <code>Step8</code>
+                                // namespace in the last line above, the
+                                // following is the main function of the
+                                // program and is again exactly like in
+                                // step-6 (apart from the changed class
+                                // names, of course).
 int main ()
 {
   try
     {
-      deallog.depth_console (0);
+      dealii::deallog.depth_console (0);
 
-      ElasticProblem<2> elastic_problem_2d;
+      Step8::ElasticProblem<2> elastic_problem_2d;
       elastic_problem_2d.run ();
     }
   catch (std::exception &exc)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.