]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move things local to each program into a local namespace.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Sep 2011 04:02:31 +0000 (04:02 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Sep 2011 04:02:31 +0000 (04:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@24295 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-40/step-40.cc
deal.II/examples/step-45/step-45.cc
deal.II/examples/step-46/step-46.cc
deal.II/examples/step-47/step-47.cc

index fcc14db60c6e4c9a15be753726f72bf90d2788b3..732097b0d9ceaccd6b4e572821c761fb13e6b14d 100644 (file)
@@ -4,7 +4,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2009, 2010 by Timo Heister and the deal.II authors */
+/*    Copyright (C) 2009, 2010, 2011 by Timo Heister and the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
 #include <fstream>
 #include <iostream>
 
-using namespace dealii;
-
-                                 // @sect3{The <code>LaplaceProblem</code> class template}
-
-                                // Next let's declare the main class of this
-                                // program. Its structure is almost exactly
-                                // that of the step-6 tutorial program. The
-                                // only significant differences are:
-                                // - The <code>mpi_communicator</code>
-                                //   variable that describes the set of
-                                //   processors we want this code to run
-                                //   on. In practice, this will be
-                                //   MPI_COMM_WORLD, i.e. all processors the
-                                //   batch scheduling system has assigned to
-                                //   this particular job.
-                                // - The presence of the <code>pcout</code>
-                                //   variable of type ConditionOStream.
-                                // - The obvious use of
-                                //   parallel::distributed::Triangulation
-                                //   instead of Triangulation.
-                                // - The presence of two IndexSet objects
-                                //   that denote which sets of degrees of
-                                //   freedom (and associated elements of
-                                //   solution and right hand side vectors) we
-                                //   own on the current processor and which
-                                //   we need (as ghost elements) for the
-                                //   algorithms in this program to work.
-                                // - The fact that all matrices and
-                                //   vectors are now distributed. We
-                                //   use their PETScWrapper versions
-                                //   for this since deal.II's own
-                                //   classes do not provide %parallel
-                                //   functionality. Note that as part
-                                //   of this class, we store a
-                                //   solution vector that does not
-                                //   only contain the degrees of
-                                //   freedom the current processor
-                                //   owns, but also (as ghost
-                                //   elements) all those vector
-                                //   elements that correspond to
-                                //   "locally relevant" degrees of
-                                //   freedom (i.e. all those that
-                                //   live on locally owned cells or
-                                //   the layer of ghost cells that
-                                //   surround it).
-template <int dim>
-class LaplaceProblem
+namespace Step40
 {
-  public:
-    LaplaceProblem ();
-    ~LaplaceProblem ();
-
-    void run ();
-
-  private:
-    void setup_system ();
-    void assemble_system ();
-    void solve ();
-    void refine_grid ();
-    void output_results (const unsigned int cycle) const;
-
-    MPI_Comm mpi_communicator;
-
-    parallel::distributed::Triangulation<dim>   triangulation;
-
-    DoFHandler<dim>      dof_handler;
-    FE_Q<dim>            fe;
-
-    IndexSet             locally_owned_dofs;
-    IndexSet             locally_relevant_dofs;
-
-    ConstraintMatrix     constraints;
+  using namespace dealii;
+
+                                  // @sect3{The <code>LaplaceProblem</code> class template}
+
+                                  // Next let's declare the main class of this
+                                  // program. Its structure is almost exactly
+                                  // that of the step-6 tutorial program. The
+                                  // only significant differences are:
+                                  // - The <code>mpi_communicator</code>
+                                  //   variable that describes the set of
+                                  //   processors we want this code to run
+                                  //   on. In practice, this will be
+                                  //   MPI_COMM_WORLD, i.e. all processors the
+                                  //   batch scheduling system has assigned to
+                                  //   this particular job.
+                                  // - The presence of the <code>pcout</code>
+                                  //   variable of type ConditionOStream.
+                                  // - The obvious use of
+                                  //   parallel::distributed::Triangulation
+                                  //   instead of Triangulation.
+                                  // - The presence of two IndexSet objects
+                                  //   that denote which sets of degrees of
+                                  //   freedom (and associated elements of
+                                  //   solution and right hand side vectors) we
+                                  //   own on the current processor and which
+                                  //   we need (as ghost elements) for the
+                                  //   algorithms in this program to work.
+                                  // - The fact that all matrices and
+                                  //   vectors are now distributed. We
+                                  //   use their PETScWrapper versions
+                                  //   for this since deal.II's own
+                                  //   classes do not provide %parallel
+                                  //   functionality. Note that as part
+                                  //   of this class, we store a
+                                  //   solution vector that does not
+                                  //   only contain the degrees of
+                                  //   freedom the current processor
+                                  //   owns, but also (as ghost
+                                  //   elements) all those vector
+                                  //   elements that correspond to
+                                  //   "locally relevant" degrees of
+                                  //   freedom (i.e. all those that
+                                  //   live on locally owned cells or
+                                  //   the layer of ghost cells that
+                                  //   surround it).
+  template <int dim>
+  class LaplaceProblem
+  {
+    public:
+      LaplaceProblem ();
+      ~LaplaceProblem ();
+
+      void run ();
+
+    private:
+      void setup_system ();
+      void assemble_system ();
+      void solve ();
+      void refine_grid ();
+      void output_results (const unsigned int cycle) const;
+
+      MPI_Comm mpi_communicator;
+
+      parallel::distributed::Triangulation<dim>   triangulation;
+
+      DoFHandler<dim>      dof_handler;
+      FE_Q<dim>            fe;
+
+      IndexSet             locally_owned_dofs;
+      IndexSet             locally_relevant_dofs;
+
+      ConstraintMatrix     constraints;
+
+      PETScWrappers::MPI::SparseMatrix system_matrix;
+      PETScWrappers::MPI::Vector locally_relevant_solution;
+      PETScWrappers::MPI::Vector system_rhs;
+
+      ConditionalOStream                pcout;
+  };
+
+
+                                  // @sect3{The <code>LaplaceProblem</code> class implementation}
+
+                                  // @sect4{Constructors and destructors}
+
+                                  // Constructors and destructors are rather
+                                  // trivial. In addition to what we do in
+                                  // step-6, we set the set of processors we
+                                  // want to work on to all machines available
+                                  // (MPI_COMM_WORLD); ask the triangulation to
+                                  // ensure that the mesh remains smooth and
+                                  // free to refined islands, for example; and
+                                  // initialize the <code>pcout</code> variable
+                                  // to only allow processor zero to output
+                                  // anything:
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem ()
+                 :
+                 mpi_communicator (MPI_COMM_WORLD),
+                 triangulation (mpi_communicator,
+                                typename Triangulation<dim>::MeshSmoothing
+                                (Triangulation<dim>::smoothing_on_refinement |
+                                 Triangulation<dim>::smoothing_on_coarsening)),
+                 dof_handler (triangulation),
+                 fe (2),
+                 pcout (std::cout,
+                        (Utilities::System::
+                         get_this_mpi_process(mpi_communicator)
+                         == 0))
+  {}
+
+
+
+  template <int dim>
+  LaplaceProblem<dim>::~LaplaceProblem ()
+  {
+    dof_handler.clear ();
+  }
+
+
+                                  // @sect4{LaplaceProblem::setup_system}
+
+                                  // The following function is, arguably, the
+                                  // most interesting one in the entire program
+                                  // since it goes to the heart of what
+                                  // distinguishes %parallel step-40 from
+                                  // sequential step-6.
+                                  //
+                                  // At the top we do what we always do: tell
+                                  // the DoFHandler object to distribute
+                                  // degrees of freedom. Since the
+                                  // triangulation we use here is distributed,
+                                  // the DoFHandler object is smart enough to
+                                  // recognize that on each processor it can
+                                  // only distribute degrees of freedom on
+                                  // cells it owns; this is followed by an
+                                  // exchange step in which processors tell
+                                  // each other about degrees of freedom on
+                                  // ghost cell. The result is a DoFHandler
+                                  // that knows about the degrees of freedom on
+                                  // locally owned cells and ghost cells
+                                  // (i.e. cells adjacent to locally owned
+                                  // cells) but nothing about cells that are
+                                  // further away, consistent with the basic
+                                  // philosophy of distributed computing that
+                                  // no processor can know everything.
+  template <int dim>
+  void LaplaceProblem<dim>::setup_system ()
+  {
+    dof_handler.distribute_dofs (fe);
+
+                                    // The next two lines extract some
+                                    // informatino we will need later
+                                    // on, namely two index sets that
+                                    // provide information about which
+                                    // degrees of freedom are owned by
+                                    // the current processor (this
+                                    // information will be used to
+                                    // initialize solution and right
+                                    // hand side vectors, and the
+                                    // system matrix, indicating which
+                                    // elements to store on the current
+                                    // processor and which to expect to
+                                    // be stored somewhere else); and
+                                    // an index set that indicates
+                                    // which degrees of freedom are
+                                    // locally relevant (i.e. live on
+                                    // cells that the current processor
+                                    // owns or on the layer of ghost
+                                    // cells around the locally owned
+                                    // cells; we need all of these
+                                    // degrees of freedom, for example,
+                                    // to estimate the error on the
+                                    // local cells).
+    locally_owned_dofs = dof_handler.locally_owned_dofs ();
+    DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                            locally_relevant_dofs);
+
+                                    // Next, let us initialize the
+                                    // solution and right hand side
+                                    // vectors. As mentioned above, the
+                                    // solution vector we seek does not
+                                    // only store elements we own, but
+                                    // also ghost entries; on the other
+                                    // hand, the right hand side vector
+                                    // only needs to have the entries
+                                    // the current processor owns since
+                                    // all we will ever do is write
+                                    // into it, never read from it on
+                                    // locally owned cells (of course
+                                    // the linear solvers will read
+                                    // from it, but they do not care
+                                    // about the geometric location of
+                                    // degrees of freedom).
+    locally_relevant_solution.reinit (mpi_communicator,
+                                     locally_owned_dofs,
+                                     locally_relevant_dofs);
+    locally_relevant_solution = 0;
+    system_rhs.reinit (mpi_communicator,
+                      dof_handler.n_dofs(),
+                      dof_handler.n_locally_owned_dofs());
+    system_rhs = 0;
+
+                                    // The next step is to compute hanging node
+                                    // and boundary value constraints, which we
+                                    // combine into a single object storing all
+                                    // constraints.
+                                    //
+                                    // As with all other things in %parallel,
+                                    // the mantra must be that no processor can
+                                    // store all information about the entire
+                                    // universe. As a consequence, we need to
+                                    // tell the constraints object for which
+                                    // degrees of freedom it can store
+                                    // constraints and for which it may not
+                                    // expect any information to store. In our
+                                    // case, as explained in the @ref
+                                    // distributed module, the degrees of
+                                    // freedom we need to care about on each
+                                    // processor are the locally relevant ones,
+                                    // so we pass this to the
+                                    // ConstraintMatrix::reinit function. As a
+                                    // side note, if you forget to pass this
+                                    // argument, the ConstraintMatrix class
+                                    // will allocate an array with length equal
+                                    // to the largest DoF index it has seen so
+                                    // far. For processors with high MPI
+                                    // process number, this may be very large
+                                    // -- maybe on the order of billions. The
+                                    // program would then allocate more memory
+                                    // than for likely all other operations
+                                    // combined for this single array.
+    constraints.clear ();
+    constraints.reinit (locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             ZeroFunction<dim>(),
+                                             constraints);
+    constraints.close ();
+
+                                    // The last part of this function deals
+                                    // with initializing the matrix with
+                                    // accompanying sparsity pattern. As in
+                                    // previous tutorial programs, we use the
+                                    // CompressedSimpleSparsityPattern as an
+                                    // intermediate with which we then
+                                    // initialize the PETSc matrix. To do so we
+                                    // have to tell the sparsity pattern its
+                                    // size but as above there is no way the
+                                    // resulting object will be able to store
+                                    // even a single pointer for each global
+                                    // degree of freedom; the best we can hope
+                                    // for is that it stores information about
+                                    // each locally relevant degree of freedom,
+                                    // i.e. all those that we may ever touch in
+                                    // the process of assembling the matrix
+                                    // (the @ref distributed_paper
+                                    // "distributed computing paper" has a long
+                                    // discussion why one really needs the
+                                    // locally relevant, and not the small set
+                                    // of locally active degrees of freedom in
+                                    // this context).
+                                    //
+                                    // So we tell the sparsity pattern its size
+                                    // and what DoFs to store anything for and
+                                    // then ask DoFTools::make_sparsity_pattern
+                                    // to fill it (this function ignores all
+                                    // cells that are not locally owned,
+                                    // mimicking what we will do below in the
+                                    // assembly process). After this, we call a
+                                    // function that exchanges entries in these
+                                    // sparsity pattern between processors so
+                                    // that in the end each processor really
+                                    // knows about all the entries that will
+                                    // exist in that part of the finite element
+                                    // matrix that it will own. The final step
+                                    // is to initialize the matrix with the
+                                    // sparsity pattern.
+    CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+                                        dof_handler.n_dofs(),
+                                        locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern (dof_handler,
+                                    csp,
+                                    constraints, false);
+    SparsityTools::distribute_sparsity_pattern (csp,
+                                               dof_handler.n_locally_owned_dofs_per_processor(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+    system_matrix.reinit (mpi_communicator,
+                         csp,
+                         dof_handler.n_locally_owned_dofs_per_processor(),
+                         dof_handler.n_locally_owned_dofs_per_processor(),
+                         Utilities::System::get_this_mpi_process(mpi_communicator));
+  }
+
+
+
+                                  // @sect4{LaplaceProblem::assemble_system}
+
+                                  // The function that then assembles the
+                                  // linear system is comparatively boring,
+                                  // being almost exactly what we've seen
+                                  // before. The points to watch out for are:
+                                  // - Assembly must only loop over locally
+                                  //   owned cells. We test this by comparing
+                                  //   a cell's subdomain_id against
+                                  //   information from the triangulation
+                                  //   but an equally valid condition would
+                                  //   have been to skip all cells for which
+                                  //   the condition <code>cell->is_ghost()
+                                  //   || cell->is_artificial()</code> is
+                                  //   true.
+                                  // - Copying local contributions into the
+                                  //   global matrix must include distributing
+                                  //   constraints and boundary values. In
+                                  //   other words, we can now (as we did in
+                                  //   step-6) first copy every local
+                                  //   contribution into the global matrix and
+                                  //   only in a later step take care of
+                                  //   hanging node constraints and boundary
+                                  //   values. The reason is, as discussed in
+                                  //   step-17, that PETSc does not provide
+                                  //   access to arbitrary elements of the
+                                  //   matrix once they have been assembled
+                                  //   into it -- in parts because they may
+                                  //   simple no longer reside on the current
+                                  //   processor but have instead been shipped
+                                  //   to a different machine.
+                                  // - The way we compute the right hand side
+                                  //   (given the formula stated in the
+                                  //   introduction) may not be the most
+                                  //   elegant but will do for a program whose
+                                  //   focus lies somewhere entirely different.
+  template <int dim>
+  void LaplaceProblem<dim>::assemble_system ()
+  {
+    const QGauss<dim>  quadrature_formula(3);
+
+    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);
+
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
+       {
+         cell_matrix = 0;
+         cell_rhs = 0;
+
+         fe_values.reinit (cell);
+
+         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+           {
+             const double
+               rhs_value
+               = (fe_values.quadrature_point(q_point)[1]
+                  >
+                  0.5+0.25*std::sin(4.0 * numbers::PI *
+                                    fe_values.quadrature_point(q_point)[0])
+                  ? 1 : -1);
+
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               {
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+                                        fe_values.shape_grad(j,q_point) *
+                                        fe_values.JxW(q_point));
+
+                 cell_rhs(i) += (rhs_value *
+                                 fe_values.shape_value(i,q_point) *
+                                 fe_values.JxW(q_point));
+               }
+           }
+
+         cell->get_dof_indices (local_dof_indices);
+         constraints.distribute_local_to_global (cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+       }
 
-    PETScWrappers::MPI::SparseMatrix system_matrix;
-    PETScWrappers::MPI::Vector locally_relevant_solution;
-    PETScWrappers::MPI::Vector system_rhs;
+    system_matrix.compress ();
+    system_rhs.compress ();
+  }
+
+
+
+                                  // @sect4{LaplaceProblem::solve}
+
+                                  // Even though solving linear systems
+                                  // on potentially tens of thousands
+                                  // of processors is by far not a
+                                  // trivial job, the function that
+                                  // does this is -- at least at the
+                                  // outside -- relatively simple. Most
+                                  // of the parts you've seen
+                                  // before. There are really only two
+                                  // things worth mentioning:
+                                  // - Solvers and preconditioners are
+                                  //   built on the deal.II wrappers of
+                                  //   PETSc functionality. It is
+                                  //   relatively well known that the
+                                  //   primary bottleneck of massively
+                                  //   %parallel linear solvers is not
+                                  //   actually the communication
+                                  //   between processors, but the fact
+                                  //   that it is difficult to produce
+                                  //   preconditioners that scale well
+                                  //   to large numbers of
+                                  //   processors. Over the second half
+                                  //   of the first decade of the 21st
+                                  //   century, it has become clear
+                                  //   that algebraic multigrid (AMG)
+                                  //   methods turn out to be extremely
+                                  //   efficient in this context, and
+                                  //   we will use one of them -- the
+                                  //   BoomerAMG implementation of the
+                                  //   Hypre package that can be
+                                  //   interfaced to through PETSc --
+                                  //   for the current program. The
+                                  //   rest of the solver itself is
+                                  //   boilerplate and has been shown
+                                  //   before. Since the linear system
+                                  //   is symmetric and positive
+                                  //   definite, we can use the CG
+                                  //   method as the outer solver.
+                                  // - Ultimately, we want a vector
+                                  //   that stores not only the
+                                  //   elements of the solution for
+                                  //   degrees of freedom the current
+                                  //   processor owns, but also all
+                                  //   other locally relevant degrees
+                                  //   of freedom. On the other hand,
+                                  //   the solver itself needs a vector
+                                  //   that is uniquely split between
+                                  //   processors, without any
+                                  //   overlap. We therefore create a
+                                  //   vector at the beginning of this
+                                  //   function that has these
+                                  //   properties, use it to solve the
+                                  //   linear system, and only assign
+                                  //   it to the vector we want at the
+                                  //   very end. This last step ensures
+                                  //   that all ghost elements are also
+                                  //   copied as necessary.
+  template <int dim>
+  void LaplaceProblem<dim>::solve ()
+  {
+    PETScWrappers::MPI::Vector
+      completely_distributed_solution (mpi_communicator,
+                                      dof_handler.n_dofs(),
+                                      dof_handler.n_locally_owned_dofs());
 
-    ConditionalOStream                pcout;
-};
+    SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
 
+    PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
 
-                                 // @sect3{The <code>LaplaceProblem</code> class implementation}
+                                    // Ask for a symmetric preconditioner by
+                                    // setting the first parameter in
+                                    // AdditionalData to true.
+    PETScWrappers::PreconditionBoomerAMG
+      preconditioner(system_matrix,
+                    PETScWrappers::PreconditionBoomerAMG::AdditionalData(true));
 
-                                 // @sect4{Constructors and destructors}
+    solver.solve (system_matrix, completely_distributed_solution, system_rhs,
+                 preconditioner);
 
-                                // Constructors and destructors are rather
-                                // trivial. In addition to what we do in
-                                // step-6, we set the set of processors we
-                                // want to work on to all machines available
-                                // (MPI_COMM_WORLD); ask the triangulation to
-                                // ensure that the mesh remains smooth and
-                                // free to refined islands, for example; and
-                                // initialize the <code>pcout</code> variable
-                                // to only allow processor zero to output
-                                // anything:
-template <int dim>
-LaplaceProblem<dim>::LaplaceProblem ()
-               :
-               mpi_communicator (MPI_COMM_WORLD),
-               triangulation (mpi_communicator,
-                              typename Triangulation<dim>::MeshSmoothing
-                              (Triangulation<dim>::smoothing_on_refinement |
-                               Triangulation<dim>::smoothing_on_coarsening)),
-               dof_handler (triangulation),
-                fe (2),
-               pcout (std::cout,
-                      (Utilities::System::
-                       get_this_mpi_process(mpi_communicator)
-                       == 0))
-{}
+    pcout << "   Solved in " << solver_control.last_step()
+         << " iterations." << std::endl;
 
+    constraints.distribute (completely_distributed_solution);
 
+    locally_relevant_solution = completely_distributed_solution;
+    locally_relevant_solution.update_ghost_values();
+  }
 
-template <int dim>
-LaplaceProblem<dim>::~LaplaceProblem ()
-{
-  dof_handler.clear ();
-}
 
 
-                                 // @sect4{LaplaceProblem::setup_system}
+                                  // @sect4{LaplaceProblem::refine_grid}
 
-                                // The following function is, arguably, the
-                                // most interesting one in the entire program
-                                // since it goes to the heart of what
-                                // distinguishes %parallel step-40 from
-                                // sequential step-6.
-                                //
-                                // At the top we do what we always do: tell
-                                // the DoFHandler object to distribute
-                                // degrees of freedom. Since the
-                                // triangulation we use here is distributed,
-                                // the DoFHandler object is smart enough to
-                                // recognize that on each processor it can
-                                // only distribute degrees of freedom on
-                                // cells it owns; this is followed by an
-                                // exchange step in which processors tell
-                                // each other about degrees of freedom on
-                                // ghost cell. The result is a DoFHandler
-                                // that knows about the degrees of freedom on
-                                // locally owned cells and ghost cells
-                                // (i.e. cells adjacent to locally owned
-                                // cells) but nothing about cells that are
-                                // further away, consistent with the basic
-                                // philosophy of distributed computing that
-                                // no processor can know everything.
-template <int dim>
-void LaplaceProblem<dim>::setup_system ()
-{
-  dof_handler.distribute_dofs (fe);
-
-                                  // The next two lines extract some
-                                  // informatino we will need later
-                                  // on, namely two index sets that
-                                  // provide information about which
-                                  // degrees of freedom are owned by
-                                  // the current processor (this
-                                  // information will be used to
-                                  // initialize solution and right
-                                  // hand side vectors, and the
-                                  // system matrix, indicating which
-                                  // elements to store on the current
-                                  // processor and which to expect to
-                                  // be stored somewhere else); and
-                                  // an index set that indicates
-                                  // which degrees of freedom are
-                                  // locally relevant (i.e. live on
-                                  // cells that the current processor
-                                  // owns or on the layer of ghost
-                                  // cells around the locally owned
-                                  // cells; we need all of these
-                                  // degrees of freedom, for example,
-                                  // to estimate the error on the
-                                  // local cells).
-  locally_owned_dofs = dof_handler.locally_owned_dofs ();
-  DoFTools::extract_locally_relevant_dofs (dof_handler,
-                                          locally_relevant_dofs);
-
-                                  // Next, let us initialize the
-                                  // solution and right hand side
-                                  // vectors. As mentioned above, the
-                                  // solution vector we seek does not
-                                  // only store elements we own, but
-                                  // also ghost entries; on the other
-                                  // hand, the right hand side vector
-                                  // only needs to have the entries
-                                  // the current processor owns since
-                                  // all we will ever do is write
-                                  // into it, never read from it on
-                                  // locally owned cells (of course
-                                  // the linear solvers will read
-                                  // from it, but they do not care
-                                  // about the geometric location of
-                                  // degrees of freedom).
-  locally_relevant_solution.reinit (mpi_communicator,
-                                   locally_owned_dofs,
-                                   locally_relevant_dofs);
-  locally_relevant_solution = 0;
-  system_rhs.reinit (mpi_communicator,
-                    dof_handler.n_dofs(),
-                    dof_handler.n_locally_owned_dofs());
-  system_rhs = 0;
-
-                                  // The next step is to compute hanging node
-                                  // and boundary value constraints, which we
-                                  // combine into a single object storing all
-                                  // constraints.
+                                  // The function that estimates the
+                                  // error and refines the grid is
+                                  // again almost exactly like the one
+                                  // in step-6. The only difference is
+                                  // that the function that flags cells
+                                  // to be refined is now in namespace
+                                  // parallel::distributed::GridRefinement
+                                  // -- a namespace that has functions
+                                  // that can communicate between all
+                                  // involved processors and determine
+                                  // global thresholds to use in
+                                  // deciding which cells to refine and
+                                  // which to coarsen.
                                   //
-                                  // As with all other things in %parallel,
-                                  // the mantra must be that no processor can
-                                  // store all information about the entire
-                                  // universe. As a consequence, we need to
-                                  // tell the constraints object for which
-                                  // degrees of freedom it can store
-                                  // constraints and for which it may not
-                                  // expect any information to store. In our
-                                  // case, as explained in the @ref
-                                  // distributed module, the degrees of
-                                  // freedom we need to care about on each
-                                  // processor are the locally relevant ones,
-                                  // so we pass this to the
-                                  // ConstraintMatrix::reinit function. As a
-                                  // side note, if you forget to pass this
-                                  // argument, the ConstraintMatrix class
-                                  // will allocate an array with length equal
-                                  // to the largest DoF index it has seen so
-                                  // far. For processors with high MPI
-                                  // process number, this may be very large
-                                  // -- maybe on the order of billions. The
-                                  // program would then allocate more memory
-                                  // than for likely all other operations
-                                  // combined for this single array.
-  constraints.clear ();
-  constraints.reinit (locally_relevant_dofs);
-  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(),
-                                           constraints);
-  constraints.close ();
-
-                                  // The last part of this function deals
-                                  // with initializing the matrix with
-                                  // accompanying sparsity pattern. As in
-                                  // previous tutorial programs, we use the
-                                  // CompressedSimpleSparsityPattern as an
-                                  // intermediate with which we then
-                                  // initialize the PETSc matrix. To do so we
-                                  // have to tell the sparsity pattern its
-                                  // size but as above there is no way the
-                                  // resulting object will be able to store
-                                  // even a single pointer for each global
-                                  // degree of freedom; the best we can hope
-                                  // for is that it stores information about
-                                  // each locally relevant degree of freedom,
-                                  // i.e. all those that we may ever touch in
-                                  // the process of assembling the matrix
-                                  // (the @ref distributed_paper  
-                                  // "distributed computing paper" has a long
-                                  // discussion why one really needs the
-                                  // locally relevant, and not the small set
-                                  // of locally active degrees of freedom in
-                                  // this context).
+                                  // Note that we didn't have to do
+                                  // anything special about the
+                                  // KellyErrorEstimator class: we just
+                                  // give it a vector with as many
+                                  // elements as the local
+                                  // triangulation has cells (locally
+                                  // owned cells, ghost cells, and
+                                  // artificial ones), but it only
+                                  // fills those entries that
+                                  // correspond to cells that are
+                                  // locally owned.
+  template <int dim>
+  void LaplaceProblem<dim>::refine_grid ()
+  {
+    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                       QGauss<dim-1>(3),
+                                       typename FunctionMap<dim>::type(),
+                                       locally_relevant_solution,
+                                       estimated_error_per_cell);
+    parallel::distributed::GridRefinement::
+      refine_and_coarsen_fixed_number (triangulation,
+                                      estimated_error_per_cell,
+                                      0.3, 0.03);
+    triangulation.execute_coarsening_and_refinement ();
+  }
+
+
+
+                                  // @sect4{LaplaceProblem::output_results}
+
+                                  // Compared to the corresponding
+                                  // function in step-6, the one here
+                                  // is a tad more complicated. There
+                                  // are two reasons: the first one is
+                                  // that we do not just want to output
+                                  // the solution but also for each
+                                  // cell which processor owns it
+                                  // (i.e. which "subdomain" it is
+                                  // in). Secondly, as discussed at
+                                  // length in step-17 and step-18,
+                                  // generating graphical data can be a
+                                  // bottleneck in parallelizing. In
+                                  // step-18, we have moved this step
+                                  // out of the actual computation but
+                                  // shifted it into a separate program
+                                  // that later combined the output
+                                  // from various processors into a
+                                  // single file. But this doesn't
+                                  // scale: if the number of processors
+                                  // is large, this may mean that the
+                                  // step of combining data on a single
+                                  // processor later becomes the
+                                  // longest running part of the
+                                  // program, or it may produce a file
+                                  // that's so large that it can't be
+                                  // visualized any more. We here
+                                  // follow a more sensible approach,
+                                  // namely creating individual files
+                                  // for each MPI process and leaving
+                                  // it to the visualization program to
+                                  // make sense of that.
                                   //
-                                  // So we tell the sparsity pattern its size
-                                  // and what DoFs to store anything for and
-                                  // then ask DoFTools::make_sparsity_pattern
-                                  // to fill it (this function ignores all
-                                  // cells that are not locally owned,
-                                  // mimicking what we will do below in the
-                                  // assembly process). After this, we call a
-                                  // function that exchanges entries in these
-                                  // sparsity pattern between processors so
-                                  // that in the end each processor really
-                                  // knows about all the entries that will
-                                  // exist in that part of the finite element
-                                  // matrix that it will own. The final step
-                                  // is to initialize the matrix with the
-                                  // sparsity pattern.
-  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
-                                      dof_handler.n_dofs(),
-                                      locally_relevant_dofs);
-  DoFTools::make_sparsity_pattern (dof_handler,
-                                  csp,
-                                  constraints, false);
-  SparsityTools::distribute_sparsity_pattern (csp,
-                                             dof_handler.n_locally_owned_dofs_per_processor(),
-                                             mpi_communicator,
-                                             locally_relevant_dofs);
-  system_matrix.reinit (mpi_communicator,
-                       csp,
-                       dof_handler.n_locally_owned_dofs_per_processor(),
-                       dof_handler.n_locally_owned_dofs_per_processor(),
-                       Utilities::System::get_this_mpi_process(mpi_communicator));
-}
-
-
-
-                                 // @sect4{LaplaceProblem::assemble_system}
-
-                                // The function that then assembles the
-                                // linear system is comparatively boring,
-                                // being almost exactly what we've seen
-                                // before. The points to watch out for are:
-                                // - Assembly must only loop over locally
-                                //   owned cells. We test this by comparing
-                                //   a cell's subdomain_id against
-                                //   information from the triangulation
-                                //   but an equally valid condition would
-                                //   have been to skip all cells for which
-                                //   the condition <code>cell->is_ghost()
-                                //   || cell->is_artificial()</code> is
-                                //   true.
-                                // - Copying local contributions into the
-                                //   global matrix must include distributing
-                                //   constraints and boundary values. In
-                                //   other words, we can now (as we did in
-                                //   step-6) first copy every local
-                                //   contribution into the global matrix and
-                                //   only in a later step take care of
-                                //   hanging node constraints and boundary
-                                //   values. The reason is, as discussed in
-                                //   step-17, that PETSc does not provide
-                                //   access to arbitrary elements of the
-                                //   matrix once they have been assembled
-                                //   into it -- in parts because they may
-                                //   simple no longer reside on the current
-                                //   processor but have instead been shipped
-                                //   to a different machine.
-                                // - The way we compute the right hand side
-                                //   (given the formula stated in the
-                                //   introduction) may not be the most
-                                //   elegant but will do for a program whose
-                                //   focus lies somewhere entirely different.
-template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
-  const QGauss<dim>  quadrature_formula(3);
-
-  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);
-
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
+                                  // To start, the top of the function
+                                  // looks like always. In addition to
+                                  // attaching the solution vector (the
+                                  // one that has entries for all
+                                  // locally relevant, not only the
+                                  // locally owned, elements), we
+                                  // attach a data vector that stores,
+                                  // for each cell, the subdomain the
+                                  // cell belongs to. This is slightly
+                                  // tricky, because of course not
+                                  // every processor knows about every
+                                  // cell. The vector we attach
+                                  // therefore has an entry for every
+                                  // cell that the current processor
+                                  // has in its mesh (locally owned
+                                  // onces, ghost cells, and artificial
+                                  // cells), but the DataOut class will
+                                  // ignore all entries that correspond
+                                  // to cells that are not owned by the
+                                  // current processor. As a
+                                  // consequence, it doesn't actually
+                                  // matter what values we write into
+                                  // these vector entries: we simply
+                                  // fill the entire vector with the
+                                  // number of the current MPI process
+                                  // (i.e. the subdomain_id of the
+                                  // current process); this correctly
+                                  // sets the values we care for,
+                                  // i.e. the entries that correspond
+                                  // to locally owned cells, while
+                                  // providing the wrong value for all
+                                  // other elements -- but these are
+                                  // then ignored anyway.
+  template <int dim>
+  void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+  {
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (locally_relevant_solution, "u");
+
+    Vector<float> subdomain (triangulation.n_active_cells());
+    for (unsigned int i=0; i<subdomain.size(); ++i)
+      subdomain(i) = triangulation.locally_owned_subdomain();
+    data_out.add_data_vector (subdomain, "subdomain");
+
+    data_out.build_patches ();
+
+                                    // The next step is to write this
+                                    // data to disk. We choose file
+                                    // names of the form
+                                    // <code>solution-XX-PPPP.vtu</code>
+                                    // where <code>XX</code> indicates
+                                    // the refinement cycle,
+                                    // <code>PPPP</code> refers to the
+                                    // processor number (enough for up
+                                    // to 10,000 processors, though we
+                                    // hope that nobody ever tries to
+                                    // generate this much data -- you
+                                    // would likely overflow all file
+                                    // system quotas), and
+                                    // <code>.vtu</code> indicates the
+                                    // XML-based Visualization Toolkit
+                                    // (VTK) file format.
+    const std::string filename = ("solution-" +
+                                 Utilities::int_to_string (cycle, 2) +
+                                 "." +
+                                 Utilities::int_to_string
+                                 (triangulation.locally_owned_subdomain(), 4));
+    std::ofstream output ((filename + ".vtu").c_str());
+    data_out.write_vtu (output);
+
+                                    // The last step is to write a
+                                    // "master record" that lists for
+                                    // the visualization program the
+                                    // names of the various files that
+                                    // combined represents the
+                                    // graphical data for the entire
+                                    // domain. The
+                                    // DataOutBase::write_pvtu_record
+                                    // does this, and it needs a list
+                                    // of filenames that we create
+                                    // first. Note that only one
+                                    // processor needs to generate this
+                                    // file; we arbitrarily choose
+                                    // processor zero to take over this
+                                    // job.
+    if (Utilities::System::get_this_mpi_process(mpi_communicator) == 0)
       {
-       cell_matrix = 0;
-       cell_rhs = 0;
-
-       fe_values.reinit (cell);
+       std::vector<std::string> filenames;
+       for (unsigned int i=0;
+            i<Utilities::System::get_n_mpi_processes(mpi_communicator);
+            ++i)
+         filenames.push_back ("solution-" +
+                              Utilities::int_to_string (cycle, 2) +
+                              "." +
+                              Utilities::int_to_string (i, 4) +
+                              ".vtu");
+
+       std::ofstream master_output ((filename + ".pvtu").c_str());
+       data_out.write_pvtu_record (master_output, filenames);
+      }
+  }
+
+
+
+                                  // @sect4{LaplaceProblem::run}
+
+                                  // The function that controls the
+                                  // overall behavior of the program is
+                                  // again like the one in step-6. The
+                                  // minor difference are the use of
+                                  // <code>pcout</code> instead of
+                                  // <code>std::cout</code> for output
+                                  // to the console (see also step-17)
+                                  // and that we only generate
+                                  // graphical output if at most 32
+                                  // processors are involved. Without
+                                  // this limit, it would be just too
+                                  // easy for people carelessly running
+                                  // this program without reading it
+                                  // first to bring down the cluster
+                                  // interconnect and fill any file
+                                  // system available :-)
+                                  //
+                                  // A functional difference to step-6
+                                  // is the use of a square domain and
+                                  // that we start with a slightly
+                                  // finer mesh (5 global refinement
+                                  // cycles) -- there just isn't much
+                                  // of a point showing a massively
+                                  // %parallel program starting on 4
+                                  // cells (although admittedly the
+                                  // point is only slightly stronger
+                                  // starting on 1024).
+  template <int dim>
+  void LaplaceProblem<dim>::run ()
+  {
+    const unsigned int n_cycles = 8;
+    for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
+      {
+       pcout << "Cycle " << cycle << ':' << std::endl;
 
-       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       if (cycle == 0)
          {
-           const double
-             rhs_value
-             = (fe_values.quadrature_point(q_point)[1]
-                >
-                0.5+0.25*std::sin(4.0 * numbers::PI *
-                                  fe_values.quadrature_point(q_point)[0])
-                ? 1 : -1);
-
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             {
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
-                 cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
-                                      fe_values.shape_grad(j,q_point) *
-                                      fe_values.JxW(q_point));
-
-               cell_rhs(i) += (rhs_value *
-                               fe_values.shape_value(i,q_point) *
-                               fe_values.JxW(q_point));
-             }
+           GridGenerator::hyper_cube (triangulation);
+           triangulation.refine_global (5);
          }
+       else
+         refine_grid ();
 
-       cell->get_dof_indices (local_dof_indices);
-       constraints.distribute_local_to_global (cell_matrix,
-                                               cell_rhs,
-                                               local_dof_indices,
-                                               system_matrix,
-                                               system_rhs);
-      }
-
-  system_matrix.compress ();
-  system_rhs.compress ();
-}
-
-
-
-                                 // @sect4{LaplaceProblem::solve}
-
-                                // Even though solving linear systems
-                                // on potentially tens of thousands
-                                // of processors is by far not a
-                                // trivial job, the function that
-                                // does this is -- at least at the
-                                // outside -- relatively simple. Most
-                                // of the parts you've seen
-                                // before. There are really only two
-                                // things worth mentioning:
-                                // - Solvers and preconditioners are
-                                //   built on the deal.II wrappers of
-                                //   PETSc functionality. It is
-                                //   relatively well known that the
-                                //   primary bottleneck of massively
-                                //   %parallel linear solvers is not
-                                //   actually the communication
-                                //   between processors, but the fact
-                                //   that it is difficult to produce
-                                //   preconditioners that scale well
-                                //   to large numbers of
-                                //   processors. Over the second half
-                                //   of the first decade of the 21st
-                                //   century, it has become clear
-                                //   that algebraic multigrid (AMG)
-                                //   methods turn out to be extremely
-                                //   efficient in this context, and
-                                //   we will use one of them -- the
-                                //   BoomerAMG implementation of the
-                                //   Hypre package that can be
-                                //   interfaced to through PETSc --
-                                //   for the current program. The
-                                //   rest of the solver itself is
-                                //   boilerplate and has been shown
-                                //   before. Since the linear system
-                                //   is symmetric and positive
-                                //   definite, we can use the CG
-                                //   method as the outer solver.
-                                // - Ultimately, we want a vector
-                                //   that stores not only the
-                                //   elements of the solution for
-                                //   degrees of freedom the current
-                                //   processor owns, but also all
-                                //   other locally relevant degrees
-                                //   of freedom. On the other hand,
-                                //   the solver itself needs a vector
-                                //   that is uniquely split between
-                                //   processors, without any
-                                //   overlap. We therefore create a
-                                //   vector at the beginning of this
-                                //   function that has these
-                                //   properties, use it to solve the
-                                //   linear system, and only assign
-                                //   it to the vector we want at the
-                                //   very end. This last step ensures
-                                //   that all ghost elements are also
-                                //   copied as necessary.
-template <int dim>
-void LaplaceProblem<dim>::solve ()
-{
-  PETScWrappers::MPI::Vector
-    completely_distributed_solution (mpi_communicator,
-                                    dof_handler.n_dofs(),
-                                    dof_handler.n_locally_owned_dofs());
-
-  SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
+       setup_system ();
 
-  PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
+       pcout << "   Number of active cells:       "
+             << triangulation.n_global_active_cells()
+             << std::endl
+             << "   Number of degrees of freedom: "
+             << dof_handler.n_dofs()
+             << std::endl;
 
-                                  // Ask for a symmetric preconditioner by
-                                  // setting the first parameter in
-                                  // AdditionalData to true.
-  PETScWrappers::PreconditionBoomerAMG
-    preconditioner(system_matrix,
-                  PETScWrappers::PreconditionBoomerAMG::AdditionalData(true));
+       assemble_system ();
+       solve ();
 
-  solver.solve (system_matrix, completely_distributed_solution, system_rhs,
-               preconditioner);
+       if (Utilities::System::get_n_mpi_processes(mpi_communicator) <= 32)
+         output_results (cycle);
 
-  pcout << "   Solved in " << solver_control.last_step()
-       << " iterations." << std::endl;
-
-  constraints.distribute (completely_distributed_solution);
-
-  locally_relevant_solution = completely_distributed_solution;
-  locally_relevant_solution.update_ghost_values();
-}
-
-
-
-                                 // @sect4{LaplaceProblem::refine_grid}
-
-                                // The function that estimates the
-                                // error and refines the grid is
-                                // again almost exactly like the one
-                                // in step-6. The only difference is
-                                // that the function that flags cells
-                                // to be refined is now in namespace
-                                // parallel::distributed::GridRefinement
-                                // -- a namespace that has functions
-                                // that can communicate between all
-                                // involved processors and determine
-                                // global thresholds to use in
-                                // deciding which cells to refine and
-                                // which to coarsen.
-                                //
-                                // Note that we didn't have to do
-                                // anything special about the
-                                // KellyErrorEstimator class: we just
-                                // give it a vector with as many
-                                // elements as the local
-                                // triangulation has cells (locally
-                                // owned cells, ghost cells, and
-                                // artificial ones), but it only
-                                // fills those entries that
-                                // correspond to cells that are
-                                // locally owned.
-template <int dim>
-void LaplaceProblem<dim>::refine_grid ()
-{
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                     QGauss<dim-1>(3),
-                                     typename FunctionMap<dim>::type(),
-                                     locally_relevant_solution,
-                                     estimated_error_per_cell);
-  parallel::distributed::GridRefinement::
-    refine_and_coarsen_fixed_number (triangulation,
-                                    estimated_error_per_cell,
-                                    0.3, 0.03);
-  triangulation.execute_coarsening_and_refinement ();
-}
-
-
-
-                                 // @sect4{LaplaceProblem::output_results}
-
-                                // Compared to the corresponding
-                                // function in step-6, the one here
-                                // is a tad more complicated. There
-                                // are two reasons: the first one is
-                                // that we do not just want to output
-                                // the solution but also for each
-                                // cell which processor owns it
-                                // (i.e. which "subdomain" it is
-                                // in). Secondly, as discussed at
-                                // length in step-17 and step-18,
-                                // generating graphical data can be a
-                                // bottleneck in parallelizing. In
-                                // step-18, we have moved this step
-                                // out of the actual computation but
-                                // shifted it into a separate program
-                                // that later combined the output
-                                // from various processors into a
-                                // single file. But this doesn't
-                                // scale: if the number of processors
-                                // is large, this may mean that the
-                                // step of combining data on a single
-                                // processor later becomes the
-                                // longest running part of the
-                                // program, or it may produce a file
-                                // that's so large that it can't be
-                                // visualized any more. We here
-                                // follow a more sensible approach,
-                                // namely creating individual files
-                                // for each MPI process and leaving
-                                // it to the visualization program to
-                                // make sense of that.
-                                //
-                                // To start, the top of the function
-                                // looks like always. In addition to
-                                // attaching the solution vector (the
-                                // one that has entries for all
-                                // locally relevant, not only the
-                                // locally owned, elements), we
-                                // attach a data vector that stores,
-                                // for each cell, the subdomain the
-                                // cell belongs to. This is slightly
-                                // tricky, because of course not
-                                // every processor knows about every
-                                // cell. The vector we attach
-                                // therefore has an entry for every
-                                // cell that the current processor
-                                // has in its mesh (locally owned
-                                // onces, ghost cells, and artificial
-                                // cells), but the DataOut class will
-                                // ignore all entries that correspond
-                                // to cells that are not owned by the
-                                // current processor. As a
-                                // consequence, it doesn't actually
-                                // matter what values we write into
-                                // these vector entries: we simply
-                                // fill the entire vector with the
-                                // number of the current MPI process
-                                // (i.e. the subdomain_id of the
-                                // current process); this correctly
-                                // sets the values we care for,
-                                // i.e. the entries that correspond
-                                // to locally owned cells, while
-                                // providing the wrong value for all
-                                // other elements -- but these are
-                                // then ignored anyway.
-template <int dim>
-void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
-{
-  DataOut<dim> data_out;
-  data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (locally_relevant_solution, "u");
-
-  Vector<float> subdomain (triangulation.n_active_cells());
-  for (unsigned int i=0; i<subdomain.size(); ++i)
-    subdomain(i) = triangulation.locally_owned_subdomain();
-  data_out.add_data_vector (subdomain, "subdomain");
-
-  data_out.build_patches ();
-
-                                  // The next step is to write this
-                                  // data to disk. We choose file
-                                  // names of the form
-                                  // <code>solution-XX-PPPP.vtu</code>
-                                  // where <code>XX</code> indicates
-                                  // the refinement cycle,
-                                  // <code>PPPP</code> refers to the
-                                  // processor number (enough for up
-                                  // to 10,000 processors, though we
-                                  // hope that nobody ever tries to
-                                  // generate this much data -- you
-                                  // would likely overflow all file
-                                  // system quotas), and
-                                  // <code>.vtu</code> indicates the
-                                  // XML-based Visualization Toolkit
-                                  // (VTK) file format.
-  const std::string filename = ("solution-" +
-                               Utilities::int_to_string (cycle, 2) +
-                               "." +
-                               Utilities::int_to_string
-                               (triangulation.locally_owned_subdomain(), 4));
-  std::ofstream output ((filename + ".vtu").c_str());
-  data_out.write_vtu (output);
-
-                                  // The last step is to write a
-                                  // "master record" that lists for
-                                  // the visualization program the
-                                  // names of the various files that
-                                  // combined represents the
-                                  // graphical data for the entire
-                                  // domain. The
-                                  // DataOutBase::write_pvtu_record
-                                  // does this, and it needs a list
-                                  // of filenames that we create
-                                  // first. Note that only one
-                                  // processor needs to generate this
-                                  // file; we arbitrarily choose
-                                  // processor zero to take over this
-                                  // job.
-  if (Utilities::System::get_this_mpi_process(mpi_communicator) == 0)
-    {
-      std::vector<std::string> filenames;
-      for (unsigned int i=0;
-          i<Utilities::System::get_n_mpi_processes(mpi_communicator);
-          ++i)
-       filenames.push_back ("solution-" +
-                            Utilities::int_to_string (cycle, 2) +
-                            "." +
-                            Utilities::int_to_string (i, 4) +
-                            ".vtu");
-
-      std::ofstream master_output ((filename + ".pvtu").c_str());
-      data_out.write_pvtu_record (master_output, filenames);
-    }
-}
-
-
-
-                                 // @sect4{LaplaceProblem::run}
-
-                                // The function that controls the
-                                // overall behavior of the program is
-                                // again like the one in step-6. The
-                                // minor difference are the use of
-                                // <code>pcout</code> instead of
-                                // <code>std::cout</code> for output
-                                // to the console (see also step-17)
-                                // and that we only generate
-                                // graphical output if at most 32
-                                // processors are involved. Without
-                                // this limit, it would be just too
-                                // easy for people carelessly running
-                                // this program without reading it
-                                // first to bring down the cluster
-                                // interconnect and fill any file
-                                // system available :-)
-                                //
-                                // A functional difference to step-6
-                                // is the use of a square domain and
-                                // that we start with a slightly
-                                // finer mesh (5 global refinement
-                                // cycles) -- there just isn't much
-                                // of a point showing a massively
-                                // %parallel program starting on 4
-                                // cells (although admittedly the
-                                // point is only slightly stronger
-                                // starting on 1024).
-template <int dim>
-void LaplaceProblem<dim>::run ()
-{
-  const unsigned int n_cycles = 8;
-  for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
-    {
-      pcout << "Cycle " << cycle << ':' << std::endl;
-
-      if (cycle == 0)
-       {
-         GridGenerator::hyper_cube (triangulation);
-         triangulation.refine_global (5);
-       }
-      else
-       refine_grid ();
-
-      setup_system ();
-
-      pcout << "   Number of active cells:       "
-           << triangulation.n_global_active_cells()
-           << std::endl
-           << "   Number of degrees of freedom: "
-           << dof_handler.n_dofs()
-           << std::endl;
-
-      assemble_system ();
-      solve ();
-
-      if (Utilities::System::get_n_mpi_processes(mpi_communicator) <= 32)
-       output_results (cycle);
-
-      pcout << std::endl;
-    }
+       pcout << std::endl;
+      }
+  }
 }
 
 
@@ -897,14 +900,17 @@ int main(int argc, char *argv[])
 {
   try
     {
+      using namespace dealii;
+      using namespace Step40;
+
       PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
       deallog.depth_console (0);
-      
+
       {
        LaplaceProblem<2> laplace_problem_2d;
        laplace_problem_2d.run ();
       }
-      
+
       PetscFinalize();
     }
   catch (std::exception &exc)
index 6a7a6b086d7c1c4493883bbbff554d6d189064f4..56494a2423c1d10e084cfe60461fb71e8bb7accf 100644 (file)
@@ -3,7 +3,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2010 by the deal.II authors */
+/*    Copyright (C) 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 <fstream>
 
 
-using namespace dealii;
-
-                                 // @sect3{The <code>LaplaceProblem</code> class}
-
-                                // The class <code>LaplaceProblem</code> is
-                                // the main class of this problem. As
-                                // mentioned in the introduction, it is
-                                // fashioned after the corresponding class in
-                                // step-3. Correspondingly, the documentation
-                                // from that tutorial program applies here as
-                                // well. The only new member variable is the
-                                // <code>constraints</code> variables that
-                                // will hold the constraints from the
-                                // periodic boundary condition. We will
-                                // initialize it in the
-                                // <code>make_periodicity_constraints()</code>
-                                // function which we call from
-                                // <code>make_grid_and_dofs()</code>.
-class LaplaceProblem
+namespace Step45
 {
-  public:
-    LaplaceProblem ();
-    void run ();
-
-  private:
-    Triangulation<2> triangulation;
-
-    FE_Q<2> fe;
-    DoFHandler<2> dof_handler;
-
-    ConstraintMatrix constraints;
-
-    SparsityPattern sparsity_pattern;
-    SparseMatrix<double> system_matrix;
-    Vector<double> system_rhs;
-    Vector<double> solution;
-
-    void assemble_system ();
-    void output_results ();
-    void make_grid_and_dofs ();
-    void make_periodicity_constraints ();
-    void solve ();   
-};
-
-
-                                 // @sect3{The <code>RightHandSide</code> class}
-
-                                // The following implements the right hand
-                                // side function discussed in the
-                                // introduction. Its implementation is
-                                // obvious given what has been shown in
-                                // step-4 before:
-class RightHandSide: public Function<2>
-{
-  public:
-    RightHandSide ();
-
-    virtual double value (const Point<2>& p,
-                         const unsigned int component = 0) const;
-};
-
-
-RightHandSide::RightHandSide ()
-               :
-               Function<2> ()
-{}
+  using namespace dealii;
+
+                                  // @sect3{The <code>LaplaceProblem</code> class}
+
+                                  // The class <code>LaplaceProblem</code> is
+                                  // the main class of this problem. As
+                                  // mentioned in the introduction, it is
+                                  // fashioned after the corresponding class in
+                                  // step-3. Correspondingly, the documentation
+                                  // from that tutorial program applies here as
+                                  // well. The only new member variable is the
+                                  // <code>constraints</code> variables that
+                                  // will hold the constraints from the
+                                  // periodic boundary condition. We will
+                                  // initialize it in the
+                                  // <code>make_periodicity_constraints()</code>
+                                  // function which we call from
+                                  // <code>make_grid_and_dofs()</code>.
+  class LaplaceProblem
+  {
+    public:
+      LaplaceProblem ();
+      void run ();
+
+    private:
+      Triangulation<2> triangulation;
+
+      FE_Q<2> fe;
+      DoFHandler<2> dof_handler;
+
+      ConstraintMatrix constraints;
+
+      SparsityPattern sparsity_pattern;
+      SparseMatrix<double> system_matrix;
+      Vector<double> system_rhs;
+      Vector<double> solution;
+
+      void assemble_system ();
+      void output_results ();
+      void make_grid_and_dofs ();
+      void make_periodicity_constraints ();
+      void solve ();
+  };
+
+
+                                  // @sect3{The <code>RightHandSide</code> class}
+
+                                  // The following implements the right hand
+                                  // side function discussed in the
+                                  // introduction. Its implementation is
+                                  // obvious given what has been shown in
+                                  // step-4 before:
+  class RightHandSide: public Function<2>
+  {
+    public:
+      RightHandSide ();
+
+      virtual double value (const Point<2>& p,
+                           const unsigned int component = 0) const;
+  };
+
+
+  RightHandSide::RightHandSide ()
+                 :
+                 Function<2> ()
+  {}
+
+
+  double
+  RightHandSide::value (const Point<2>&p,
+                       const unsigned int) const
+  {
+    return (std::cos (2 * numbers::PI * p(0)) *
+           std::exp (- 2 * p(0)) *
+           std::cos (2 * numbers::PI * p(1)) *
+           std::exp (- 2 * p(1)));
+  }
+
+                                  // @sect3{Implementation of the <code>LaplaceProblem</code> class}
+
+                                  // The first part of implementing the main
+                                  // class is the constructor. It is unchanged
+                                  // from step-3 and step-4:
+  LaplaceProblem::LaplaceProblem ()
+                 :
+                 fe (1),
+                 dof_handler (triangulation)
+  {}
+
+
+                                  // @sect4{LaplaceProblem::make_grid_and_dofs}
+
+                                  // The following is the first function to be
+                                  // called in <code>run()</code>. It sets up
+                                  // the mesh and degrees of freedom.
+                                  //
+                                  // We start by creating the usual square mesh
+                                  // and changing the boundary indicator on the
+                                  // parts of the boundary where we have
+                                  // Dirichlet boundary conditions (top and
+                                  // bottom, i.e. faces two and three of the
+                                  // reference cell as defined by
+                                  // GeometryInfo), so that we can distinguish
+                                  // between the parts of the boundary where
+                                  // periodic and where Dirichlet boundary
+                                  // conditions hold. We then refine the mesh a
+                                  // fixed number of times, with child faces
+                                  // inheriting the boundary indicators
+                                  // previously set on the coarse mesh from
+                                  // their parents.
+  void LaplaceProblem::make_grid_and_dofs ()
+  {
+    GridGenerator::hyper_cube (triangulation);
+    triangulation.begin_active ()->face (2)->set_boundary_indicator (1);
+    triangulation.begin_active ()->face (3)->set_boundary_indicator (1);
+    triangulation.refine_global (5);
+
+                                    // The next step is to distribute the
+                                    // degrees of freedom and produce a little
+                                    // bit of graphical output:
+    dof_handler.distribute_dofs (fe);
+    std::cout << "Number of active cells: "
+             << triangulation.n_active_cells ()
+             << std::endl
+             << "Degrees of freedom: " << dof_handler.n_dofs ()
+             << std::endl;
+
+                                    // Now it is the time for the constraints
+                                    // that come from the periodicity
+                                    // constraints. We do this in the
+                                    // following, separate function, after
+                                    // clearing any possible prior content from
+                                    // the constraints object:
+    constraints.clear ();
+    make_periodicity_constraints ();
+
+                                    // We also incorporate the homogeneous
+                                    // Dirichlet boundary conditions on the
+                                    // upper and lower parts of the boundary
+                                    // (i.e. the ones with boundary indicator
+                                    // 1) and close the
+                                    // <code>ConstraintMatrix</code> object:
+    VectorTools::interpolate_boundary_values (dof_handler, 1,
+                                             ZeroFunction<2> (),
+                                             constraints);
+    constraints.close ();
+
+                                    // Then we create the sparsity pattern and
+                                    // the system matrix and initialize the
+                                    // solution and right-hand side
+                                    // vectors. This is again as in step-3 or
+                                    // step-6, for example:
+    CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs(),
+                                                 dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler,
+                                    c_sparsity_pattern,
+                                    constraints,
+                                    false);
+    c_sparsity_pattern.compress ();
+    sparsity_pattern.copy_from (c_sparsity_pattern);
+
+    system_matrix.reinit (sparsity_pattern);
+    system_rhs.reinit (dof_handler.n_dofs());
+    solution.reinit (dof_handler.n_dofs());
+  }
+
+
+
+                                  // @sect4{LaplaceProblem::make_periodicity_constraints}
+
+                                  // This is the function that provides the new
+                                  // material of this tutorial program. The
+                                  // general outline of the algorithm is as
+                                  // follows: we first loop over all the
+                                  // degrees of freedom on the right boundary
+                                  // and record their $y$-locations in a map
+                                  // together with their global indices. Then
+                                  // we go along the left boundary, find
+                                  // matching $y$-locations for each degree of
+                                  // freedom, and then add constraints that
+                                  // identify these matched degrees of freedom.
+                                  //
+                                  // In this function, we make use of the fact
+                                  // that we have a scalar element (i.e. the
+                                  // only valid vector component that can be
+                                  // passed to DoFAccessor::vertex_dof_index is
+                                  // zero) and that we have a $Q_1$ element for
+                                  // which all degrees of freedom live in the
+                                  // vertices of the cell. Furthermore, we have
+                                  // assumed that we are in 2d and that meshes
+                                  // were not refined adaptively &mdash; the
+                                  // latter assumption would imply that there
+                                  // may be vertices that aren't matched
+                                  // one-to-one and for which we won't be able
+                                  // to compute constraints this easily. We
+                                  // will discuss in the "outlook" part of the
+                                  // results section below other strategies to
+                                  // write the current function that can work
+                                  // in cases like this as well.
+  void LaplaceProblem::make_periodicity_constraints ()
+  {
+                                    // To start with the actual implementation,
+                                    // we loop over all active cells and check
+                                    // whether the cell is located at the right
+                                    // boundary (i.e. face 1 &mdash; the one at
+                                    // the right end of the cell &mdash; is at
+                                    // the boundary). If that is so, then we
+                                    // use that for the currently used finite
+                                    // element, each degree of freedom of the
+                                    // face is located on one vertex, and store
+                                    // their $y$-coordinate along with the
+                                    // global number of this degree of freedom
+                                    // in the following map:
+    std::map<unsigned int, double> dof_locations;
+
+    for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
+        cell != dof_handler.end (); ++cell)
+      if (cell->at_boundary ()
+         &&
+         cell->face(1)->at_boundary ())
+       {
+         dof_locations[cell->face(1)->vertex_dof_index(0, 0)]
+           = cell->face(1)->vertex(0)[1];
+         dof_locations[cell->face(1)->vertex_dof_index(1, 0)]
+           = cell->face(1)->vertex(1)[1];
+       }
+                                    // Note that in the above block, we add
+                                    // vertices zero and one of the affected
+                                    // face to the map. This means that we will
+                                    // add each vertex twice, once from each of
+                                    // the two adjacent cells (unless the
+                                    // vertex is a corner of the domain). Since
+                                    // the coordinates of the vertex are the
+                                    // same both times of course, there is no
+                                    // harm: we replace one value in the map
+                                    // with itself the second time we visit an
+                                    // entry.
+                                    //
+                                    // The same will be true below where we add
+                                    // the same constraint twice to the
+                                    // ConstraintMatrix &mdash; again, we will
+                                    // overwrite the constraint with itself,
+                                    // and no harm is done.
+
+                                    // Now we have to find the corresponding
+                                    // degrees of freedom on the left part of
+                                    // the boundary. Therefore we loop over all
+                                    // cells again and choose the ones where
+                                    // face 0 is at the boundary:
+    for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
+        cell != dof_handler.end (); ++cell)
+      if (cell->at_boundary ()
+         &&
+         cell->face (0)->at_boundary ())
+       {
+                                          // Every degree of freedom on this
+                                          // face needs to have a corresponding
+                                          // one on the right side of the face,
+                                          // and our goal is to add a
+                                          // constraint for the one on the left
+                                          // in terms of the one on the
+                                          // right. To this end we first add a
+                                          // new line to the constraint matrix
+                                          // for this one degree of
+                                          // freedom. Then we identify it with
+                                          // the corresponding degree of
+                                          // freedom on the right part of the
+                                          // boundary by constraining the
+                                          // degree of freedom on the left with
+                                          // the one on the right times a
+                                          // weight of 1.0.
+                                          //
+                                          // Consequently, we loop over the two
+                                          // vertices of each face we find and
+                                          // then loop over all the
+                                          // $y$-locations we've previously
+                                          // recorded to find which degree of
+                                          // freedom on the right boundary
+                                          // corresponds to the one we
+                                          // currently look at. Note that we
+                                          // have entered these into a map, and
+                                          // when looping over the iterators
+                                          // <code>p</code> of this map,
+                                          // <code>p-@>first</code> corresponds
+                                          // to the "key" of an entry (the
+                                          // global number of the degree of
+                                          // freedom), whereas
+                                          // <code>p-@>second</code> is the
+                                          // "value" (the $y$-location we have
+                                          // entered above).
+                                          //
+                                          // We are quite sure here that we
+                                          // should be finding such a
+                                          // corresponding degree of
+                                          // freedom. However, sometimes stuff
+                                          // happens and so the bottom of the
+                                          // block contains an assertion that
+                                          // our assumption was indeed correct
+                                          // and that a vertex was found.
+         for (unsigned int face_vertex = 0; face_vertex<2; ++face_vertex)
+           {
+             constraints.add_line (cell->face(0)->vertex_dof_index (face_vertex, 0));
+
+             std::map<unsigned int, double>::const_iterator p = dof_locations.begin();
+             for (; p != dof_locations.end(); ++p)
+               if (std::fabs(p->second - cell->face(0)->vertex(face_vertex)[1]) < 1e-8)
+                 {
+                   constraints.add_entry (cell->face(0)->vertex_dof_index (face_vertex, 0),
+                                          p->first, 1.0);
+                   break;
+                 }
+             Assert (p != dof_locations.end(),
+                     ExcMessage ("No corresponding degree of freedom was found!"));
+           }
+       }
+  }
+
+
+
+                                  // @sect4{LaplaceProblem::assemble_system}
+
+                                  // Assembling the system matrix and the
+                                  // right-hand side vector is done as in other
+                                  // tutorials before.
+                                  //
+                                  // The only difference here is that we don't
+                                  // copy elements from local contributions
+                                  // into the global matrix and later fix up
+                                  // constrained degrees of freedom, but that
+                                  // we let the ConstraintMatrix do this job in
+                                  // one swoop for us using the
+                                  // ConstraintMatrix::distribute_local_to_global
+                                  // function(). This was previously already
+                                  // demonstrated in step-16, step-22, for
+                                  // example, along with a discussion in the
+                                  // introduction of step-27.
+  void LaplaceProblem::assemble_system ()
+  {
+    QGauss<2>  quadrature_formula(2);
+    FEValues<2> 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();
 
-double
-RightHandSide::value (const Point<2>&p,
-                     const unsigned int) const
-{
-  return (std::cos (2 * numbers::PI * p(0)) *
-         std::exp (- 2 * p(0)) *
-         std::cos (2 * numbers::PI * p(1)) *
-         std::exp (- 2 * p(1)));
-}
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       cell_rhs (dofs_per_cell);
 
-                                 // @sect3{Implementation of the <code>LaplaceProblem</code> class}
-
-                                // The first part of implementing the main
-                                // class is the constructor. It is unchanged
-                                // from step-3 and step-4:
-LaplaceProblem::LaplaceProblem ()
-               :
-               fe (1),
-               dof_handler (triangulation)
-{}
-
-
-                                 // @sect4{LaplaceProblem::make_grid_and_dofs}
-
-                                // The following is the first function to be
-                                // called in <code>run()</code>. It sets up
-                                // the mesh and degrees of freedom.
-                                //
-                                // We start by creating the usual square mesh
-                                // and changing the boundary indicator on the
-                                // parts of the boundary where we have
-                                // Dirichlet boundary conditions (top and
-                                // bottom, i.e. faces two and three of the
-                                // reference cell as defined by
-                                // GeometryInfo), so that we can distinguish
-                                // between the parts of the boundary where
-                                // periodic and where Dirichlet boundary
-                                // conditions hold. We then refine the mesh a
-                                // fixed number of times, with child faces
-                                // inheriting the boundary indicators
-                                // previously set on the coarse mesh from
-                                // their parents.
-void LaplaceProblem::make_grid_and_dofs ()
-{
-  GridGenerator::hyper_cube (triangulation);
-  triangulation.begin_active ()->face (2)->set_boundary_indicator (1);
-  triangulation.begin_active ()->face (3)->set_boundary_indicator (1);   
-  triangulation.refine_global (5);
-  
-                                  // The next step is to distribute the
-                                  // degrees of freedom and produce a little
-                                  // bit of graphical output:
-  dof_handler.distribute_dofs (fe);
-  std::cout << "Number of active cells: "
-           << triangulation.n_active_cells ()
-           << std::endl
-           << "Degrees of freedom: " << dof_handler.n_dofs ()
-           << std::endl;
-  
-                                  // Now it is the time for the constraints
-                                  // that come from the periodicity
-                                  // constraints. We do this in the
-                                  // following, separate function, after
-                                  // clearing any possible prior content from
-                                  // the constraints object:
-  constraints.clear ();
-  make_periodicity_constraints ();
-
-                                  // We also incorporate the homogeneous
-                                  // Dirichlet boundary conditions on the
-                                  // upper and lower parts of the boundary
-                                  // (i.e. the ones with boundary indicator
-                                  // 1) and close the
-                                  // <code>ConstraintMatrix</code> object:
-  VectorTools::interpolate_boundary_values (dof_handler, 1,
-                                           ZeroFunction<2> (),
-                                           constraints);
-  constraints.close ();
-  
-                                  // Then we create the sparsity pattern and
-                                  // the system matrix and initialize the
-                                  // solution and right-hand side
-                                  // vectors. This is again as in step-3 or
-                                  // step-6, for example:
-  CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs(),
-                                               dof_handler.n_dofs());
-  DoFTools::make_sparsity_pattern (dof_handler,
-                                  c_sparsity_pattern,
-                                  constraints,
-                                  false);
-  c_sparsity_pattern.compress ();
-  sparsity_pattern.copy_from (c_sparsity_pattern);
-  
-  system_matrix.reinit (sparsity_pattern);
-  system_rhs.reinit (dof_handler.n_dofs());
-  solution.reinit (dof_handler.n_dofs());
-}
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
+    const RightHandSide right_hand_side;
 
-
-                                 // @sect4{LaplaceProblem::make_periodicity_constraints}
-
-                                // This is the function that provides the new
-                                // material of this tutorial program. The
-                                // general outline of the algorithm is as
-                                // follows: we first loop over all the
-                                // degrees of freedom on the right boundary
-                                // and record their $y$-locations in a map
-                                // together with their global indices. Then
-                                // we go along the left boundary, find
-                                // matching $y$-locations for each degree of
-                                // freedom, and then add constraints that
-                                // identify these matched degrees of freedom.
-                                //
-                                // In this function, we make use of the fact
-                                // that we have a scalar element (i.e. the
-                                // only valid vector component that can be
-                                // passed to DoFAccessor::vertex_dof_index is
-                                // zero) and that we have a $Q_1$ element for
-                                // which all degrees of freedom live in the
-                                // vertices of the cell. Furthermore, we have
-                                // assumed that we are in 2d and that meshes
-                                // were not refined adaptively &mdash; the
-                                // latter assumption would imply that there
-                                // may be vertices that aren't matched
-                                // one-to-one and for which we won't be able
-                                // to compute constraints this easily. We
-                                // will discuss in the "outlook" part of the
-                                // results section below other strategies to
-                                // write the current function that can work
-                                // in cases like this as well.
-void LaplaceProblem::make_periodicity_constraints ()
-{
-                                  // To start with the actual implementation,
-                                  // we loop over all active cells and check
-                                  // whether the cell is located at the right
-                                  // boundary (i.e. face 1 &mdash; the one at
-                                  // the right end of the cell &mdash; is at
-                                  // the boundary). If that is so, then we
-                                  // use that for the currently used finite
-                                  // element, each degree of freedom of the
-                                  // face is located on one vertex, and store
-                                  // their $y$-coordinate along with the
-                                  // global number of this degree of freedom
-                                  // in the following map:
-  std::map<unsigned int, double> dof_locations;
-
-  for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
-       cell != dof_handler.end (); ++cell)
-    if (cell->at_boundary ()
-       &&
-       cell->face(1)->at_boundary ())
-      {
-       dof_locations[cell->face(1)->vertex_dof_index(0, 0)]
-         = cell->face(1)->vertex(0)[1];
-       dof_locations[cell->face(1)->vertex_dof_index(1, 0)]
-         = cell->face(1)->vertex(1)[1];
-      }
-                                  // Note that in the above block, we add
-                                  // vertices zero and one of the affected
-                                  // face to the map. This means that we will
-                                  // add each vertex twice, once from each of
-                                  // the two adjacent cells (unless the
-                                  // vertex is a corner of the domain). Since
-                                  // the coordinates of the vertex are the
-                                  // same both times of course, there is no
-                                  // harm: we replace one value in the map
-                                  // with itself the second time we visit an
-                                  // entry.
-                                  //
-                                  // The same will be true below where we add
-                                  // the same constraint twice to the
-                                  // ConstraintMatrix &mdash; again, we will
-                                  // overwrite the constraint with itself,
-                                  // and no harm is done.
-  
-                                  // Now we have to find the corresponding
-                                  // degrees of freedom on the left part of
-                                  // the boundary. Therefore we loop over all
-                                  // cells again and choose the ones where
-                                  // face 0 is at the boundary:
-  for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
-       cell != dof_handler.end (); ++cell)
-    if (cell->at_boundary ()
-       &&
-       cell->face (0)->at_boundary ())
+    DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
+                                       endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
       {
-                                        // Every degree of freedom on this
-                                        // face needs to have a corresponding
-                                        // one on the right side of the face,
-                                        // and our goal is to add a
-                                        // constraint for the one on the left
-                                        // in terms of the one on the
-                                        // right. To this end we first add a
-                                        // new line to the constraint matrix
-                                        // for this one degree of
-                                        // freedom. Then we identify it with
-                                        // the corresponding degree of
-                                        // freedom on the right part of the
-                                        // boundary by constraining the
-                                        // degree of freedom on the left with
-                                        // the one on the right times a
-                                        // weight of 1.0.
-                                        //
-                                        // Consequently, we loop over the two
-                                        // vertices of each face we find and
-                                        // then loop over all the
-                                        // $y$-locations we've previously
-                                        // recorded to find which degree of
-                                        // freedom on the right boundary
-                                        // corresponds to the one we
-                                        // currently look at. Note that we
-                                        // have entered these into a map, and
-                                        // when looping over the iterators
-                                        // <code>p</code> of this map,
-                                        // <code>p-@>first</code> corresponds
-                                        // to the "key" of an entry (the
-                                        // global number of the degree of
-                                        // freedom), whereas
-                                        // <code>p-@>second</code> is the
-                                        // "value" (the $y$-location we have
-                                        // entered above).
-                                        //
-                                        // We are quite sure here that we
-                                        // should be finding such a
-                                        // corresponding degree of
-                                        // freedom. However, sometimes stuff
-                                        // happens and so the bottom of the
-                                        // block contains an assertion that
-                                        // our assumption was indeed correct
-                                        // and that a vertex was found.
-       for (unsigned int face_vertex = 0; face_vertex<2; ++face_vertex)
-         {
-           constraints.add_line (cell->face(0)->vertex_dof_index (face_vertex, 0));
-
-           std::map<unsigned int, double>::const_iterator p = dof_locations.begin();
-           for (; p != dof_locations.end(); ++p)
-             if (std::fabs(p->second - cell->face(0)->vertex(face_vertex)[1]) < 1e-8)
-               {
-                 constraints.add_entry (cell->face(0)->vertex_dof_index (face_vertex, 0),
-                                        p->first, 1.0);
-                 break;
-               }
-           Assert (p != dof_locations.end(),
-                   ExcMessage ("No corresponding degree of freedom was found!"));
-         }
+       fe_values.reinit (cell);
+       cell_matrix = 0;
+       cell_rhs = 0;
+
+       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           {
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
+                                    fe_values.shape_grad (j, q_point) *
+                                    fe_values.JxW (q_point));
+
+             cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+                             right_hand_side.value (fe_values.quadrature_point (q_point)) *
+                             fe_values.JxW (q_point));
+           }
+
+       cell->get_dof_indices (local_dof_indices);
+       constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                               local_dof_indices,
+                                               system_matrix, system_rhs);
       }
-}
-
-
+  }
 
-                                 // @sect4{LaplaceProblem::assemble_system}
-
-                                // Assembling the system matrix and the
-                                // right-hand side vector is done as in other
-                                // tutorials before.
-                                //
-                                // The only difference here is that we don't
-                                // copy elements from local contributions
-                                // into the global matrix and later fix up
-                                // constrained degrees of freedom, but that
-                                // we let the ConstraintMatrix do this job in
-                                // one swoop for us using the
-                                // ConstraintMatrix::distribute_local_to_global
-                                // function(). This was previously already
-                                // demonstrated in step-16, step-22, for
-                                // example, along with a discussion in the
-                                // introduction of step-27.
-void LaplaceProblem::assemble_system ()
-{
-  QGauss<2>  quadrature_formula(2);
-  FEValues<2> 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();
+                                  // @sect4{LaplaceProblem::solve}
 
-  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double>       cell_rhs (dofs_per_cell);
+                                  // To solve the linear system of equations
+                                  // $Au=b$ we use the CG solver with an
+                                  // SSOR-preconditioner. This is, again,
+                                  // copied almost verbatim from step-4, with
+                                  // the exception of the preconditioner. As in
+                                  // step-6, we need to make sure that
+                                  // constrained degrees of freedom get their
+                                  // correct values after solving by calling
+                                  // the ConstraintMatrix::distribute function:
+  void LaplaceProblem::solve ()
+  {
+    SolverControl solver_control (dof_handler.n_dofs (), 1e-12);
+    PreconditionSSOR<SparseMatrix<double> > precondition;
 
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    precondition.initialize (system_matrix);
 
-  const RightHandSide right_hand_side;
+    SolverCG<> cg (solver_control);
 
-  DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
-                                     endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      fe_values.reinit (cell);
-      cell_matrix = 0;
-      cell_rhs = 0;
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-       for (unsigned int i=0; i<dofs_per_cell; ++i)
-         {
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
-                                  fe_values.shape_grad (j, q_point) *
-                                  fe_values.JxW (q_point));
-
-           cell_rhs(i) += (fe_values.shape_value (i, q_point) *
-                           right_hand_side.value (fe_values.quadrature_point (q_point)) *
-                           fe_values.JxW (q_point));
-         }
-      
-      cell->get_dof_indices (local_dof_indices);
-      constraints.distribute_local_to_global (cell_matrix, cell_rhs,
-                                             local_dof_indices,
-                                             system_matrix, system_rhs);
-    }
-}
+    cg.solve (system_matrix, solution, system_rhs, precondition);
+    constraints.distribute (solution);
+  }
 
 
-                                 // @sect4{LaplaceProblem::solve}
+                                  // @sect4{LaplaceProblem::output_results}
 
-                                // To solve the linear system of equations
-                                // $Au=b$ we use the CG solver with an
-                                // SSOR-preconditioner. This is, again,
-                                // copied almost verbatim from step-4, with
-                                // the exception of the preconditioner. As in
-                                // step-6, we need to make sure that
-                                // constrained degrees of freedom get their
-                                // correct values after solving by calling
-                                // the ConstraintMatrix::distribute function:
-void LaplaceProblem::solve ()
-{
-  SolverControl solver_control (dof_handler.n_dofs (), 1e-12);
-  PreconditionSSOR<SparseMatrix<double> > precondition;
-   
-  precondition.initialize (system_matrix);
-   
-  SolverCG<> cg (solver_control);
-   
-  cg.solve (system_matrix, solution, system_rhs, precondition);
-  constraints.distribute (solution);
-}
+                                  // This is another function copied from
+                                  // previous tutorial programs. It generates
+                                  // graphical output in VTK format:
+  void LaplaceProblem::output_results ()
+  {
+    DataOut<2> data_out;
 
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, "u");
+    data_out.build_patches ();
 
-                                 // @sect4{LaplaceProblem::output_results}
+    std::ofstream output ("solution.vtk");
 
-                                // This is another function copied from
-                                // previous tutorial programs. It generates
-                                // graphical output in VTK format:
-void LaplaceProblem::output_results ()
-{
-  DataOut<2> data_out;
-   
-  data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (solution, "u");
-  data_out.build_patches ();
-   
-  std::ofstream output ("solution.vtk");
-   
-  data_out.write_vtk (output);
-}
+    data_out.write_vtk (output);
+  }
 
 
 
-                                 // @sect4{LaplaceProblem::run}
+                                  // @sect4{LaplaceProblem::run}
 
-                                // And another function copied from previous
-                                // programs:
-void LaplaceProblem::run ()
-{
-  make_grid_and_dofs();
-  assemble_system ();
-  solve ();
-  output_results ();
+                                  // And another function copied from previous
+                                  // programs:
+  void LaplaceProblem::run ()
+  {
+    make_grid_and_dofs();
+    assemble_system ();
+    solve ();
+    output_results ();
+  }
 }
 
                                  // @sect3{The <code>main</code> function}
@@ -491,6 +494,9 @@ int main ()
 {
   try
     {
+      using namespace dealii;
+      using namespace Step45;
+
       deallog.depth_console (0);
 
       LaplaceProblem laplace_problem;
@@ -509,7 +515,7 @@ int main ()
 
       return 1;
     }
-  catch (...) 
+  catch (...)
     {
       std::cerr << std::endl << std::endl
                << "----------------------------------------------------"
index 52a0882115d4a083efaf3381e5b76a316d56d936..75dc091cfbf6565155fbc6617e95b1668108300e 100644 (file)
 #include <fstream>
 #include <sstream>
 
-using namespace dealii;
-
-                                 // @sect3{The <code>FluidStructureProblem</code> class template}
-
-                                // This is the main class. It is, if you
-                                // want, a combination of step-8 and step-22
-                                // in that it has member variables that
-                                // either address the global problem (the
-                                // Triangulation and hp::DoFHandler objects,
-                                // as well as the hp::FECollection and
-                                // various linear algebra objects) or that
-                                // pertain to either the elasticity or Stokes
-                                // sub-problems. The general structure of the
-                                // class, however, is like that of most of
-                                // the other programs implementing stationary
-                                // problems.
-                                //
-                                // There are a few helper functions
-                                // (<code>cell_is_in_fluid_domain,
-                                // cell_is_in_solid_domain</code>) of
-                                // self-explanatory nature (operating on the
-                                // symbolic names for the two subdomains that
-                                // will be used as material_ids for cells
-                                // belonging to the subdomains, as explained
-                                // in the introduction) and a few functions
-                                // (<code>make_grid, set_active_fe_indices,
-                                // assemble_interface_terms</code>) that have
-                                // been broken out of other functions that
-                                // can be found in many of the other tutorial
-                                // programs and that will be discussed as we
-                                // get to their implementation.
-                                //
-                                // The final set of variables
-                                // (<code>viscosity, lambda, eta</code>)
-                                // describes the material properties used for
-                                // the two physics models.
-template <int dim>
-class FluidStructureProblem
-{
-  public:
-    FluidStructureProblem (const unsigned int stokes_degree,
-                          const unsigned int elasticity_degree);
-    void run ();
-
-  private:
-    enum
-    {
-         fluid_domain_id,
-         solid_domain_id
-    };
-
-    static bool
-    cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
-
-    static bool
-    cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
-
-
-    void make_grid ();
-    void set_active_fe_indices ();
-    void setup_dofs ();
-    void assemble_system ();
-    void assemble_interface_term (const FEFaceValuesBase<dim>          &elasticity_fe_face_values,
-                                 const FEFaceValuesBase<dim>          &stokes_fe_face_values,
-                                 std::vector<Tensor<1,dim> >          &elasticity_phi,
-                                 std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
-                                 std::vector<double>                  &stokes_phi_p,
-                                 FullMatrix<double>                   &local_interface_matrix) const;
-    void solve ();
-    void output_results (const unsigned int refinement_cycle) const;
-    void refine_mesh ();
-
-    const unsigned int    stokes_degree;
-    const unsigned int    elasticity_degree;
-
-    Triangulation<dim>    triangulation;
-    FESystem<dim>         stokes_fe;
-    FESystem<dim>         elasticity_fe;
-    hp::FECollection<dim> fe_collection;
-    hp::DoFHandler<dim>   dof_handler;
-
-    ConstraintMatrix      constraints;
-
-    SparsityPattern       sparsity_pattern;
-    SparseMatrix<double>  system_matrix;
-
-    Vector<double>        solution;
-    Vector<double>        system_rhs;
-
-    const double          viscosity;
-    const double          lambda;
-    const double          mu;
-};
-
-
-                                 // @sect3{Boundary values and right hand side}
-
-                                // The following classes do as their names
-                                // suggest. The boundary values for the
-                                // velocity are $\mathbf u=(0, \sin(\pi
-                                // x))^T$ in 2d and $\mathbf u=(0, 0,
-                                // \sin(\pi x)\sin(\pi y))^T$ in 3d,
-                                // respectively. The remaining boundary
-                                // conditions for this problem are all
-                                // homogenous and have been discussed in the
-                                // introduction. The right hand side forcing
-                                // term is zero for both the fluid and the
-                                // solid.
-template <int dim>
-class StokesBoundaryValues : public Function<dim>
-{
-  public:
-    StokesBoundaryValues () : Function<dim>(dim+1+dim) {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
-};
-
 
-template <int dim>
-double
-StokesBoundaryValues<dim>::value (const Point<dim>  &p,
-                                 const unsigned int component) const
+namespace Step46
 {
-  Assert (component < this->n_components,
-         ExcIndexRange (component, 0, this->n_components));
+  using namespace dealii;
+
+                                  // @sect3{The <code>FluidStructureProblem</code> class template}
+
+                                  // This is the main class. It is, if you
+                                  // want, a combination of step-8 and step-22
+                                  // in that it has member variables that
+                                  // either address the global problem (the
+                                  // Triangulation and hp::DoFHandler objects,
+                                  // as well as the hp::FECollection and
+                                  // various linear algebra objects) or that
+                                  // pertain to either the elasticity or Stokes
+                                  // sub-problems. The general structure of the
+                                  // class, however, is like that of most of
+                                  // the other programs implementing stationary
+                                  // problems.
+                                  //
+                                  // There are a few helper functions
+                                  // (<code>cell_is_in_fluid_domain,
+                                  // cell_is_in_solid_domain</code>) of
+                                  // self-explanatory nature (operating on the
+                                  // symbolic names for the two subdomains that
+                                  // will be used as material_ids for cells
+                                  // belonging to the subdomains, as explained
+                                  // in the introduction) and a few functions
+                                  // (<code>make_grid, set_active_fe_indices,
+                                  // assemble_interface_terms</code>) that have
+                                  // been broken out of other functions that
+                                  // can be found in many of the other tutorial
+                                  // programs and that will be discussed as we
+                                  // get to their implementation.
+                                  //
+                                  // The final set of variables
+                                  // (<code>viscosity, lambda, eta</code>)
+                                  // describes the material properties used for
+                                  // the two physics models.
+  template <int dim>
+  class FluidStructureProblem
+  {
+    public:
+      FluidStructureProblem (const unsigned int stokes_degree,
+                            const unsigned int elasticity_degree);
+      void run ();
 
-  if (component == dim-1)
-    switch (dim)
+    private:
+      enum
       {
-       case 2:
-             return std::sin(numbers::PI*p[0]);
-       case 3:
-             return std::sin(numbers::PI*p[0]) * std::sin(numbers::PI*p[1]);
-       default:
-             Assert (false, ExcNotImplemented());
-      }
+           fluid_domain_id,
+           solid_domain_id
+      };
+
+      static bool
+      cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
+
+      static bool
+      cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
+
+
+      void make_grid ();
+      void set_active_fe_indices ();
+      void setup_dofs ();
+      void assemble_system ();
+      void assemble_interface_term (const FEFaceValuesBase<dim>          &elasticity_fe_face_values,
+                                   const FEFaceValuesBase<dim>          &stokes_fe_face_values,
+                                   std::vector<Tensor<1,dim> >          &elasticity_phi,
+                                   std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
+                                   std::vector<double>                  &stokes_phi_p,
+                                   FullMatrix<double>                   &local_interface_matrix) const;
+      void solve ();
+      void output_results (const unsigned int refinement_cycle) const;
+      void refine_mesh ();
+
+      const unsigned int    stokes_degree;
+      const unsigned int    elasticity_degree;
+
+      Triangulation<dim>    triangulation;
+      FESystem<dim>         stokes_fe;
+      FESystem<dim>         elasticity_fe;
+      hp::FECollection<dim> fe_collection;
+      hp::DoFHandler<dim>   dof_handler;
+
+      ConstraintMatrix      constraints;
+
+      SparsityPattern       sparsity_pattern;
+      SparseMatrix<double>  system_matrix;
+
+      Vector<double>        solution;
+      Vector<double>        system_rhs;
+
+      const double          viscosity;
+      const double          lambda;
+      const double          mu;
+  };
+
+
+                                  // @sect3{Boundary values and right hand side}
+
+                                  // The following classes do as their names
+                                  // suggest. The boundary values for the
+                                  // velocity are $\mathbf u=(0, \sin(\pi
+                                  // x))^T$ in 2d and $\mathbf u=(0, 0,
+                                  // \sin(\pi x)\sin(\pi y))^T$ in 3d,
+                                  // respectively. The remaining boundary
+                                  // conditions for this problem are all
+                                  // homogenous and have been discussed in the
+                                  // introduction. The right hand side forcing
+                                  // term is zero for both the fluid and the
+                                  // solid.
+  template <int dim>
+  class StokesBoundaryValues : public Function<dim>
+  {
+    public:
+      StokesBoundaryValues () : Function<dim>(dim+1+dim) {}
 
-  return 0;
-}
+      virtual double value (const Point<dim>   &p,
+                           const unsigned int  component = 0) const;
 
+      virtual void vector_value (const Point<dim> &p,
+                                Vector<double>   &value) const;
+  };
 
-template <int dim>
-void
-StokesBoundaryValues<dim>::vector_value (const Point<dim> &p,
-                                        Vector<double>   &values) const
-{
-  for (unsigned int c=0; c<this->n_components; ++c)
-    values(c) = StokesBoundaryValues<dim>::value (p, c);
-}
 
+  template <int dim>
+  double
+  StokesBoundaryValues<dim>::value (const Point<dim>  &p,
+                                   const unsigned int component) const
+  {
+    Assert (component < this->n_components,
+           ExcIndexRange (component, 0, this->n_components));
 
+    if (component == dim-1)
+      switch (dim)
+       {
+         case 2:
+               return std::sin(numbers::PI*p[0]);
+         case 3:
+               return std::sin(numbers::PI*p[0]) * std::sin(numbers::PI*p[1]);
+         default:
+               Assert (false, ExcNotImplemented());
+       }
 
-template <int dim>
-class RightHandSide : public Function<dim>
-{
-  public:
-    RightHandSide () : Function<dim>(dim+1) {}
+    return 0;
+  }
 
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
 
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
+  template <int dim>
+  void
+  StokesBoundaryValues<dim>::vector_value (const Point<dim> &p,
+                                          Vector<double>   &values) const
+  {
+    for (unsigned int c=0; c<this->n_components; ++c)
+      values(c) = StokesBoundaryValues<dim>::value (p, c);
+  }
 
-};
 
 
-template <int dim>
-double
-RightHandSide<dim>::value (const Point<dim>  &/*p*/,
-                           const unsigned int /*component*/) const
-{
-  return 0;
-}
+  template <int dim>
+  class RightHandSide : public Function<dim>
+  {
+    public:
+      RightHandSide () : Function<dim>(dim+1) {}
 
+      virtual double value (const Point<dim>   &p,
+                           const unsigned int  component = 0) const;
 
-template <int dim>
-void
-RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                  Vector<double>   &values) const
-{
-  for (unsigned int c=0; c<this->n_components; ++c)
-    values(c) = RightHandSide<dim>::value (p, c);
-}
+      virtual void vector_value (const Point<dim> &p,
+                                Vector<double>   &value) const;
 
+  };
 
 
-                                 // @sect3{The <code>FluidStructureProblem</code> implementation}
-
-                                 // @sect4{Constructors and helper functions}
-
-                                // Let's now get to the implementation of the
-                                // primary class of this program. The first
-                                // few functions are the constructor and the
-                                // helper functions that can be used to
-                                // determine which part of the domain a cell
-                                // is in. Given the discussion of these
-                                // topics in the introduction, their
-                                // implementation is rather obvious. In the
-                                // constructor, note that we have to
-                                // construct the hp::FECollection object from
-                                // the base elements for Stokes and
-                                // elasticity; using the
-                                // hp::FECollection::push_back function
-                                // assigns them spots zero and one in this
-                                // collection, an order that we have to
-                                // remember and use consistently in the rest
-                                // of the program.
-template <int dim>
-FluidStructureProblem<dim>::
-FluidStructureProblem (const unsigned int stokes_degree,
-                      const unsigned int elasticity_degree)
-                :
-                stokes_degree (stokes_degree),
-               elasticity_degree (elasticity_degree),
-                triangulation (Triangulation<dim>::maximum_smoothing),
-                stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
-                          FE_Q<dim>(stokes_degree), 1,
-                          FE_Nothing<dim>(), dim),
-                elasticity_fe (FE_Nothing<dim>(), dim,
-                              FE_Nothing<dim>(), 1,
-                              FE_Q<dim>(elasticity_degree), dim),
-                dof_handler (triangulation),
-                viscosity (2),
-               lambda (1),
-               mu (1)
-{
-  fe_collection.push_back (stokes_fe);
-  fe_collection.push_back (elasticity_fe);
-}
+  template <int dim>
+  double
+  RightHandSide<dim>::value (const Point<dim>  &/*p*/,
+                            const unsigned int /*component*/) const
+  {
+    return 0;
+  }
 
 
+  template <int dim>
+  void
+  RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                   Vector<double>   &values) const
+  {
+    for (unsigned int c=0; c<this->n_components; ++c)
+      values(c) = RightHandSide<dim>::value (p, c);
+  }
 
 
-template <int dim>
-bool
-FluidStructureProblem<dim>::
-cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
-{
-  return (cell->material_id() == fluid_domain_id);
-}
 
+                                  // @sect3{The <code>FluidStructureProblem</code> implementation}
+
+                                  // @sect4{Constructors and helper functions}
+
+                                  // Let's now get to the implementation of the
+                                  // primary class of this program. The first
+                                  // few functions are the constructor and the
+                                  // helper functions that can be used to
+                                  // determine which part of the domain a cell
+                                  // is in. Given the discussion of these
+                                  // topics in the introduction, their
+                                  // implementation is rather obvious. In the
+                                  // constructor, note that we have to
+                                  // construct the hp::FECollection object from
+                                  // the base elements for Stokes and
+                                  // elasticity; using the
+                                  // hp::FECollection::push_back function
+                                  // assigns them spots zero and one in this
+                                  // collection, an order that we have to
+                                  // remember and use consistently in the rest
+                                  // of the program.
+  template <int dim>
+  FluidStructureProblem<dim>::
+  FluidStructureProblem (const unsigned int stokes_degree,
+                        const unsigned int elasticity_degree)
+                 :
+                 stokes_degree (stokes_degree),
+                 elasticity_degree (elasticity_degree),
+                 triangulation (Triangulation<dim>::maximum_smoothing),
+                 stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
+                            FE_Q<dim>(stokes_degree), 1,
+                            FE_Nothing<dim>(), dim),
+                 elasticity_fe (FE_Nothing<dim>(), dim,
+                                FE_Nothing<dim>(), 1,
+                                FE_Q<dim>(elasticity_degree), dim),
+                 dof_handler (triangulation),
+                 viscosity (2),
+                 lambda (1),
+                 mu (1)
+  {
+    fe_collection.push_back (stokes_fe);
+    fe_collection.push_back (elasticity_fe);
+  }
 
-template <int dim>
-bool
-FluidStructureProblem<dim>::
-cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
-{
-  return (cell->material_id() == solid_domain_id);
-}
 
 
-                                 // @sect4{Meshes and assigning subdomains}
-
-                                // The next pair of functions deals with
-                                // generating a mesh and making sure all
-                                // flags that denote subdomains are
-                                // correct. <code>make_grid</code>, as
-                                // discussed in the introduction, generates
-                                // an $8\times 8$ mesh (or an $8\times
-                                // 8\times 8$ mesh in 3d) to make sure that
-                                // each coarse mesh cell is completely within
-                                // one of the subdomains. After generating
-                                // this mesh, we loop over its boundary and
-                                // set the boundary indicator to one at the
-                                // top boundary, the only place where we set
-                                // nonzero Dirichlet boundary
-                                // conditions. After this, we loop again over
-                                // all cells to set the material indicator
-                                // &mdash; used to denote which part of the
-                                // domain we are in, to either the fluid or
-                                // solid indicator.
-template <int dim>
-void
-FluidStructureProblem<dim>::make_grid ()
-{
-  GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1);
-
-  for (typename Triangulation<dim>::active_cell_iterator
-        cell = triangulation.begin_active();
-       cell != triangulation.end(); ++cell)
-    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-      if (cell->face(f)->at_boundary()
-         &&
-         (cell->face(f)->center()[dim-1] == 1))
-       cell->face(f)->set_all_boundary_indicators(1);
-
-
-  for (typename Triangulation<dim>::active_cell_iterator
-         cell = dof_handler.begin_active();
-       cell != dof_handler.end(); ++cell)
-    if (((std::fabs(cell->center()[0]) < 0.25)
-         &&
-         (cell->center()[dim-1] > 0.5))
-       ||
-       ((std::fabs(cell->center()[0]) >= 0.25)
-        &&
-        (cell->center()[dim-1] > -0.5)))
-      cell->set_material_id (fluid_domain_id);
-    else
-      cell->set_material_id (solid_domain_id);
-}
 
+  template <int dim>
+  bool
+  FluidStructureProblem<dim>::
+  cell_is_in_fluid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
+  {
+    return (cell->material_id() == fluid_domain_id);
+  }
 
-                                // The second part of this pair of functions
-                                // determines which finite element to use on
-                                // each cell. Above we have set the material
-                                // indicator for each coarse mesh cell, and
-                                // as mentioned in the introduction, this
-                                // information is inherited from mother to
-                                // child cell upon mesh refinement.
-                                //
-                                // In other words, whenever we have refined
-                                // (or created) the mesh, we can rely on the
-                                // material indicators to be a correct
-                                // description of which part of the domain a
-                                // cell is in. We then use this to set the
-                                // active FE index of the cell to the
-                                // corresponding element of the
-                                // hp::FECollection member variable of this
-                                // class: zero for fluid cells, one for solid
-                                // cells.
-template <int dim>
-void
-FluidStructureProblem<dim>::set_active_fe_indices ()
-{
-  for (typename hp::DoFHandler<dim>::active_cell_iterator
-         cell = dof_handler.begin_active();
-       cell != dof_handler.end(); ++cell)
-    {
-      if (cell_is_in_fluid_domain(cell))
-       cell->set_active_fe_index (0);
-      else if (cell_is_in_solid_domain(cell))
-       cell->set_active_fe_index (1);
-      else
-       Assert (false, ExcNotImplemented());
-    }
-}
 
+  template <int dim>
+  bool
+  FluidStructureProblem<dim>::
+  cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell)
+  {
+    return (cell->material_id() == solid_domain_id);
+  }
 
-                                 // @sect4{<code>FluidStructureProblem::setup_dofs</code>}
-
-                                // The next step is to setup the data
-                                // structures for the linear system. To this
-                                // end, we first have to set the active FE
-                                // indices with the function immediately
-                                // above, then distribute degrees of freedom,
-                                // and then determine constraints on the
-                                // linear system. The latter includes hanging
-                                // node constraints as usual, but also the
-                                // inhomogenous boundary values at the top
-                                // fluid boundary, and zero boundary values
-                                // along the perimeter of the solid
-                                // subdomain.
-template <int dim>
-void
-FluidStructureProblem<dim>::setup_dofs ()
-{
-  set_active_fe_indices ();
-  dof_handler.distribute_dofs (fe_collection);
 
+                                  // @sect4{Meshes and assigning subdomains}
+
+                                  // The next pair of functions deals with
+                                  // generating a mesh and making sure all
+                                  // flags that denote subdomains are
+                                  // correct. <code>make_grid</code>, as
+                                  // discussed in the introduction, generates
+                                  // an $8\times 8$ mesh (or an $8\times
+                                  // 8\times 8$ mesh in 3d) to make sure that
+                                  // each coarse mesh cell is completely within
+                                  // one of the subdomains. After generating
+                                  // this mesh, we loop over its boundary and
+                                  // set the boundary indicator to one at the
+                                  // top boundary, the only place where we set
+                                  // nonzero Dirichlet boundary
+                                  // conditions. After this, we loop again over
+                                  // all cells to set the material indicator
+                                  // &mdash; used to denote which part of the
+                                  // domain we are in, to either the fluid or
+                                  // solid indicator.
+  template <int dim>
+  void
+  FluidStructureProblem<dim>::make_grid ()
   {
-    constraints.clear ();
-    DoFTools::make_hanging_node_constraints (dof_handler,
-                                            constraints);
+    GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1);
 
-    std::vector<bool> velocity_mask (dim+1+dim, false);
-    for (unsigned int d=0; d<dim; ++d)
-      velocity_mask[d] = true;
-    VectorTools::interpolate_boundary_values (dof_handler,
-                                             1,
-                                             StokesBoundaryValues<dim>(),
-                                             constraints,
-                                             velocity_mask);
-
-    std::vector<bool> elasticity_mask (dim+1+dim, false);
-    for (unsigned int d=dim+1; d<dim+1+dim; ++d)
-      elasticity_mask[d] = true;
-    VectorTools::interpolate_boundary_values (dof_handler,
-                                             0,
-                                             ZeroFunction<dim>(dim+1+dim),
-                                             constraints,
-                                             elasticity_mask);
+    for (typename Triangulation<dim>::active_cell_iterator
+          cell = triangulation.begin_active();
+        cell != triangulation.end(); ++cell)
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+       if (cell->face(f)->at_boundary()
+           &&
+           (cell->face(f)->center()[dim-1] == 1))
+         cell->face(f)->set_all_boundary_indicators(1);
+
+
+    for (typename Triangulation<dim>::active_cell_iterator
+          cell = dof_handler.begin_active();
+        cell != dof_handler.end(); ++cell)
+      if (((std::fabs(cell->center()[0]) < 0.25)
+          &&
+          (cell->center()[dim-1] > 0.5))
+         ||
+         ((std::fabs(cell->center()[0]) >= 0.25)
+          &&
+          (cell->center()[dim-1] > -0.5)))
+       cell->set_material_id (fluid_domain_id);
+      else
+       cell->set_material_id (solid_domain_id);
   }
 
-                                  // There are more constraints we have to
-                                  // handle, though: we have to make sure
-                                  // that the velocity is zero at the
-                                  // interface between fluid and solid. The
-                                  // following piece of code was already
-                                  // presented in the introduction:
+
+                                  // The second part of this pair of functions
+                                  // determines which finite element to use on
+                                  // each cell. Above we have set the material
+                                  // indicator for each coarse mesh cell, and
+                                  // as mentioned in the introduction, this
+                                  // information is inherited from mother to
+                                  // child cell upon mesh refinement.
+                                  //
+                                  // In other words, whenever we have refined
+                                  // (or created) the mesh, we can rely on the
+                                  // material indicators to be a correct
+                                  // description of which part of the domain a
+                                  // cell is in. We then use this to set the
+                                  // active FE index of the cell to the
+                                  // corresponding element of the
+                                  // hp::FECollection member variable of this
+                                  // class: zero for fluid cells, one for solid
+                                  // cells.
+  template <int dim>
+  void
+  FluidStructureProblem<dim>::set_active_fe_indices ()
   {
-    std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
     for (typename hp::DoFHandler<dim>::active_cell_iterator
           cell = dof_handler.begin_active();
         cell != dof_handler.end(); ++cell)
-      if (cell_is_in_fluid_domain (cell))
-       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-         if (!cell->at_boundary(f))
-           {
-             bool face_is_on_interface = false;
-
-             if ((cell->neighbor(f)->has_children() == false)
-                 &&
-                 (cell_is_in_solid_domain (cell->neighbor(f))))
-               face_is_on_interface = true;
-             else if (cell->neighbor(f)->has_children() == true)
-               {
-                 for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
-                   if (cell_is_in_solid_domain (cell->neighbor_child_on_subface
-                                                (f, sf)))
-                     {
-                       face_is_on_interface = true;
-                       break;
-                     }
-               }
-
-             if (face_is_on_interface)
-               {
-                 cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
-                 for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
-                   if (stokes_fe.face_system_to_component_index(i).first < dim)
-                     constraints.add_line (local_face_dof_indices[i]);
-               }
-           }
+      {
+       if (cell_is_in_fluid_domain(cell))
+         cell->set_active_fe_index (0);
+       else if (cell_is_in_solid_domain(cell))
+         cell->set_active_fe_index (1);
+       else
+         Assert (false, ExcNotImplemented());
+      }
   }
 
-                                  // At the end of all this, we can declare
-                                  // to the constraints object that we now
-                                  // have all constraints ready to go and
-                                  // that the object can rebuild its internal
-                                  // data structures for better efficiency:
-  constraints.close ();
-
-  std::cout << "   Number of active cells: "
-            << triangulation.n_active_cells()
-            << std::endl
-            << "   Number of degrees of freedom: "
-            << dof_handler.n_dofs()
-            << std::endl;
-
-                                  // In the rest of this function we create a
-                                  // sparsity pattern as discussed
-                                  // extensively in the introduction, and use
-                                  // it to initialize the matrix; then also
-                                  // set vectors to their correct sizes:
+
+                                  // @sect4{<code>FluidStructureProblem::setup_dofs</code>}
+
+                                  // The next step is to setup the data
+                                  // structures for the linear system. To this
+                                  // end, we first have to set the active FE
+                                  // indices with the function immediately
+                                  // above, then distribute degrees of freedom,
+                                  // and then determine constraints on the
+                                  // linear system. The latter includes hanging
+                                  // node constraints as usual, but also the
+                                  // inhomogenous boundary values at the top
+                                  // fluid boundary, and zero boundary values
+                                  // along the perimeter of the solid
+                                  // subdomain.
+  template <int dim>
+  void
+  FluidStructureProblem<dim>::setup_dofs ()
   {
-    CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
-                                        dof_handler.n_dofs());
+    set_active_fe_indices ();
+    dof_handler.distribute_dofs (fe_collection);
+
+    {
+      constraints.clear ();
+      DoFTools::make_hanging_node_constraints (dof_handler,
+                                              constraints);
+
+      std::vector<bool> velocity_mask (dim+1+dim, false);
+      for (unsigned int d=0; d<dim; ++d)
+       velocity_mask[d] = true;
+      VectorTools::interpolate_boundary_values (dof_handler,
+                                               1,
+                                               StokesBoundaryValues<dim>(),
+                                               constraints,
+                                               velocity_mask);
+
+      std::vector<bool> elasticity_mask (dim+1+dim, false);
+      for (unsigned int d=dim+1; d<dim+1+dim; ++d)
+       elasticity_mask[d] = true;
+      VectorTools::interpolate_boundary_values (dof_handler,
+                                               0,
+                                               ZeroFunction<dim>(dim+1+dim),
+                                               constraints,
+                                               elasticity_mask);
+    }
 
-    Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(),
-                                              fe_collection.n_components());
-    Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(),
-                                              fe_collection.n_components());
+                                    // There are more constraints we have to
+                                    // handle, though: we have to make sure
+                                    // that the velocity is zero at the
+                                    // interface between fluid and solid. The
+                                    // following piece of code was already
+                                    // presented in the introduction:
+    {
+      std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
+      for (typename hp::DoFHandler<dim>::active_cell_iterator
+            cell = dof_handler.begin_active();
+          cell != dof_handler.end(); ++cell)
+       if (cell_is_in_fluid_domain (cell))
+         for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+           if (!cell->at_boundary(f))
+             {
+               bool face_is_on_interface = false;
+
+               if ((cell->neighbor(f)->has_children() == false)
+                   &&
+                   (cell_is_in_solid_domain (cell->neighbor(f))))
+                 face_is_on_interface = true;
+               else if (cell->neighbor(f)->has_children() == true)
+                 {
+                   for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
+                     if (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+                                                  (f, sf)))
+                       {
+                         face_is_on_interface = true;
+                         break;
+                       }
+                 }
 
-    for (unsigned int c=0; c<fe_collection.n_components(); ++c)
-      for (unsigned int d=0; d<fe_collection.n_components(); ++d)
-       {
-         if (((c<dim+1) && (d<dim+1)
-              && !((c==dim) && (d==dim)))
-             ||
-             ((c>=dim+1) && (d>=dim+1)))
-           cell_coupling[c][d] = DoFTools::always;
-
-         if ((c>=dim+1) && (d<dim+1))
-           face_coupling[c][d] = DoFTools::always;
-       }
-    
-    DoFTools::make_flux_sparsity_pattern (dof_handler, csp,
-                                         cell_coupling, face_coupling);
-    constraints.condense (csp);
-    sparsity_pattern.copy_from (csp);
+               if (face_is_on_interface)
+                 {
+                   cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
+                   for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
+                     if (stokes_fe.face_system_to_component_index(i).first < dim)
+                       constraints.add_line (local_face_dof_indices[i]);
+                 }
+             }
+    }
+
+                                    // At the end of all this, we can declare
+                                    // to the constraints object that we now
+                                    // have all constraints ready to go and
+                                    // that the object can rebuild its internal
+                                    // data structures for better efficiency:
+    constraints.close ();
+
+    std::cout << "   Number of active cells: "
+             << triangulation.n_active_cells()
+             << std::endl
+             << "   Number of degrees of freedom: "
+             << dof_handler.n_dofs()
+             << std::endl;
+
+                                    // In the rest of this function we create a
+                                    // sparsity pattern as discussed
+                                    // extensively in the introduction, and use
+                                    // it to initialize the matrix; then also
+                                    // set vectors to their correct sizes:
+    {
+      CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+                                          dof_handler.n_dofs());
+
+      Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(),
+                                                fe_collection.n_components());
+      Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(),
+                                                fe_collection.n_components());
+
+      for (unsigned int c=0; c<fe_collection.n_components(); ++c)
+       for (unsigned int d=0; d<fe_collection.n_components(); ++d)
+         {
+           if (((c<dim+1) && (d<dim+1)
+                && !((c==dim) && (d==dim)))
+               ||
+               ((c>=dim+1) && (d>=dim+1)))
+             cell_coupling[c][d] = DoFTools::always;
+
+           if ((c>=dim+1) && (d<dim+1))
+             face_coupling[c][d] = DoFTools::always;
+         }
+
+      DoFTools::make_flux_sparsity_pattern (dof_handler, csp,
+                                           cell_coupling, face_coupling);
+      constraints.condense (csp);
+      sparsity_pattern.copy_from (csp);
+    }
+
+    system_matrix.reinit (sparsity_pattern);
+
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
   }
 
-  system_matrix.reinit (sparsity_pattern);
 
-  solution.reinit (dof_handler.n_dofs());
-  system_rhs.reinit (dof_handler.n_dofs());
-}
 
+                                  // @sect4{<code>FluidStructureProblem::assemble_system</code>}
 
+                                  // Following is the central function of this
+                                  // program: the one that assembles the linear
+                                  // system. It has a long section of setting
+                                  // up auxiliary functions at the beginning:
+                                  // from creating the quadrature formulas and
+                                  // setting up the FEValues, FEFaceValues and
+                                  // FESubfaceValues objects necessary to
+                                  // integrate the cell terms as well as the
+                                  // interface terms for the case where cells
+                                  // along the interface come together at same
+                                  // size or with differing levels of
+                                  // refinement...
+  template <int dim>
+  void FluidStructureProblem<dim>::assemble_system ()
+  {
+    system_matrix=0;
+    system_rhs=0;
+
+    const QGauss<dim> stokes_quadrature(stokes_degree+2);
+    const QGauss<dim> elasticity_quadrature(elasticity_degree+2);
+
+    hp::QCollection<dim>  q_collection;
+    q_collection.push_back (stokes_quadrature);
+    q_collection.push_back (elasticity_quadrature);
+
+    hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
+                                   update_values    |
+                                   update_quadrature_points  |
+                                   update_JxW_values |
+                                   update_gradients);
+
+    const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
+                                                        elasticity_degree+2));
+
+    FEFaceValues<dim>    stokes_fe_face_values (stokes_fe,
+                                               common_face_quadrature,
+                                               update_JxW_values |
+                                               update_normal_vectors |
+                                               update_gradients);
+    FEFaceValues<dim>    elasticity_fe_face_values (elasticity_fe,
+                                                   common_face_quadrature,
+                                                   update_values);
+    FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
+                                                  common_face_quadrature,
+                                                  update_JxW_values |
+                                                  update_normal_vectors |
+                                                  update_gradients);
+    FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe,
+                                                      common_face_quadrature,
+                                                      update_values);
+
+                                    // ...to objects that are needed to
+                                    // describe the local contributions to the
+                                    // global linear system...
+    const unsigned int        stokes_dofs_per_cell     = stokes_fe.dofs_per_cell;
+    const unsigned int        elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
+
+    FullMatrix<double>        local_matrix;
+    FullMatrix<double>        local_interface_matrix (elasticity_dofs_per_cell,
+                                                     stokes_dofs_per_cell);
+    Vector<double>            local_rhs;
+
+    std::vector<unsigned int> local_dof_indices;
+    std::vector<unsigned int> neighbor_dof_indices (stokes_dofs_per_cell);
+
+    const RightHandSide<dim>  right_hand_side;
+
+                                    // ...to variables that allow us to extract
+                                    // certain components of the shape
+                                    // functions and cache their values rather
+                                    // than having to recompute them at every
+                                    // quadrature point:
+    const FEValuesExtractors::Vector     velocities (0);
+    const FEValuesExtractors::Scalar     pressure (dim);
+    const FEValuesExtractors::Vector     displacements (dim+1);
+
+    std::vector<SymmetricTensor<2,dim> > stokes_phi_grads_u (stokes_dofs_per_cell);
+    std::vector<double>                  stokes_div_phi_u   (stokes_dofs_per_cell);
+    std::vector<double>                  stokes_phi_p       (stokes_dofs_per_cell);
+
+    std::vector<Tensor<2,dim> >          elasticity_phi_grad (elasticity_dofs_per_cell);
+    std::vector<double>                  elasticity_phi_div  (elasticity_dofs_per_cell);
+    std::vector<Tensor<1,dim> >          elasticity_phi      (elasticity_dofs_per_cell);
+
+                                    // Then comes the main loop over all cells
+                                    // and, as in step-27, the initialization
+                                    // of the hp::FEValues object for the
+                                    // current cell and the extraction of a
+                                    // FEValues object that is appropriate for
+                                    // the current cell:
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+       hp_fe_values.reinit (cell);
+
+       const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+
+       local_matrix.reinit (cell->get_fe().dofs_per_cell,
+                            cell->get_fe().dofs_per_cell);
+       local_rhs.reinit (cell->get_fe().dofs_per_cell);
+
+                                        // With all of this done, we continue
+                                        // to assemble the cell terms for cells
+                                        // that are part of the Stokes and
+                                        // elastic regions. While we could in
+                                        // principle do this in one formula, in
+                                        // effect implementing the one bilinear
+                                        // form stated in the introduction, we
+                                        // realize that our finite element
+                                        // spaces are chosen in such a way that
+                                        // on each cell, one set of variables
+                                        // (either velocities and pressure, or
+                                        // displacements) are always zero, and
+                                        // consequently a more efficient way of
+                                        // computing local integrals is to do
+                                        // only what's necessary based on an
+                                        // <code>if</code> clause that tests
+                                        // which part of the domain we are in.
+                                        //
+                                        // The actual computation of the local
+                                        // matrix is the same as in step-22 as
+                                        // well as that given in the @ref
+                                        // vector_valued documentation module
+                                        // for the elasticity equations:
+       if (cell_is_in_fluid_domain (cell))
+         {
+           const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+           Assert (dofs_per_cell == stokes_dofs_per_cell,
+                   ExcInternalError());
 
-                                 // @sect4{<code>FluidStructureProblem::assemble_system</code>}
+           for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+             {
+               for (unsigned int k=0; k<dofs_per_cell; ++k)
+                 {
+                   stokes_phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+                   stokes_div_phi_u[k]   = fe_values[velocities].divergence (k, q);
+                   stokes_phi_p[k]       = fe_values[pressure].value (k, q);
+                 }
 
-                                // Following is the central function of this
-                                // program: the one that assembles the linear
-                                // system. It has a long section of setting
-                                // up auxiliary functions at the beginning:
-                                // from creating the quadrature formulas and
-                                // setting up the FEValues, FEFaceValues and
-                                // FESubfaceValues objects necessary to
-                                // integrate the cell terms as well as the
-                                // interface terms for the case where cells
-                                // along the interface come together at same
-                                // size or with differing levels of
-                                // refinement...
-template <int dim>
-void FluidStructureProblem<dim>::assemble_system ()
-{
-  system_matrix=0;
-  system_rhs=0;
-
-  const QGauss<dim> stokes_quadrature(stokes_degree+2);
-  const QGauss<dim> elasticity_quadrature(elasticity_degree+2);
-
-  hp::QCollection<dim>  q_collection;
-  q_collection.push_back (stokes_quadrature);
-  q_collection.push_back (elasticity_quadrature);
-
-  hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
-                                 update_values    |
-                                 update_quadrature_points  |
-                                 update_JxW_values |
-                                 update_gradients);
-
-  const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
-                                                      elasticity_degree+2));
-
-  FEFaceValues<dim>    stokes_fe_face_values (stokes_fe,
-                                             common_face_quadrature,
-                                             update_JxW_values |
-                                             update_normal_vectors |
-                                             update_gradients);
-  FEFaceValues<dim>    elasticity_fe_face_values (elasticity_fe,
-                                                 common_face_quadrature,
-                                                 update_values);
-  FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
-                                                common_face_quadrature,
-                                                update_JxW_values |
-                                                update_normal_vectors |
-                                                update_gradients);
-  FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe,
-                                                    common_face_quadrature,
-                                                    update_values);
-
-                                  // ...to objects that are needed to
-                                  // describe the local contributions to the
-                                  // global linear system...
-  const unsigned int        stokes_dofs_per_cell     = stokes_fe.dofs_per_cell;
-  const unsigned int        elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
-
-  FullMatrix<double>        local_matrix;
-  FullMatrix<double>        local_interface_matrix (elasticity_dofs_per_cell,
-                                                   stokes_dofs_per_cell);
-  Vector<double>            local_rhs;
-
-  std::vector<unsigned int> local_dof_indices;
-  std::vector<unsigned int> neighbor_dof_indices (stokes_dofs_per_cell);
-
-  const RightHandSide<dim>  right_hand_side;
-
-                                  // ...to variables that allow us to extract
-                                  // certain components of the shape
-                                  // functions and cache their values rather
-                                  // than having to recompute them at every
-                                  // quadrature point:
-  const FEValuesExtractors::Vector     velocities (0);
-  const FEValuesExtractors::Scalar     pressure (dim);
-  const FEValuesExtractors::Vector     displacements (dim+1);
-
-  std::vector<SymmetricTensor<2,dim> > stokes_phi_grads_u (stokes_dofs_per_cell);
-  std::vector<double>                  stokes_div_phi_u   (stokes_dofs_per_cell);
-  std::vector<double>                  stokes_phi_p       (stokes_dofs_per_cell);
-
-  std::vector<Tensor<2,dim> >          elasticity_phi_grad (elasticity_dofs_per_cell);
-  std::vector<double>                  elasticity_phi_div  (elasticity_dofs_per_cell);
-  std::vector<Tensor<1,dim> >          elasticity_phi      (elasticity_dofs_per_cell);
-
-                                  // Then comes the main loop over all cells
-                                  // and, as in step-27, the initialization
-                                  // of the hp::FEValues object for the
-                                  // current cell and the extraction of a
-                                  // FEValues object that is appropriate for
-                                  // the current cell:
-  typename hp::DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      hp_fe_values.reinit (cell);
-
-      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
-
-      local_matrix.reinit (cell->get_fe().dofs_per_cell,
-                          cell->get_fe().dofs_per_cell);
-      local_rhs.reinit (cell->get_fe().dofs_per_cell);
-
-                                      // With all of this done, we continue
-                                      // to assemble the cell terms for cells
-                                      // that are part of the Stokes and
-                                      // elastic regions. While we could in
-                                      // principle do this in one formula, in
-                                      // effect implementing the one bilinear
-                                      // form stated in the introduction, we
-                                      // realize that our finite element
-                                      // spaces are chosen in such a way that
-                                      // on each cell, one set of variables
-                                      // (either velocities and pressure, or
-                                      // displacements) are always zero, and
-                                      // consequently a more efficient way of
-                                      // computing local integrals is to do
-                                      // only what's necessary based on an
-                                      // <code>if</code> clause that tests
-                                      // which part of the domain we are in.
-                                      //
-                                      // The actual computation of the local
-                                      // matrix is the same as in step-22 as
-                                      // well as that given in the @ref
-                                      // vector_valued documentation module
-                                      // for the elasticity equations:
-      if (cell_is_in_fluid_domain (cell))
-       {
-         const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
-         Assert (dofs_per_cell == stokes_dofs_per_cell,
-                 ExcInternalError());
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   local_matrix(i,j) += (2 * viscosity * stokes_phi_grads_u[i] * stokes_phi_grads_u[j]
+                                         - stokes_div_phi_u[i] * stokes_phi_p[j]
+                                         - stokes_phi_p[i] * stokes_div_phi_u[j])
+                                        * fe_values.JxW(q);
+             }
+         }
+       else
+         {
+           const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+           Assert (dofs_per_cell == elasticity_dofs_per_cell,
+                   ExcInternalError());
 
-         for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
-           {
-             for (unsigned int k=0; k<dofs_per_cell; ++k)
-               {
-                 stokes_phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
-                 stokes_div_phi_u[k]   = fe_values[velocities].divergence (k, q);
-                 stokes_phi_p[k]       = fe_values[pressure].value (k, q);
-               }
-
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
-                 local_matrix(i,j) += (2 * viscosity * stokes_phi_grads_u[i] * stokes_phi_grads_u[j]
-                                       - stokes_div_phi_u[i] * stokes_phi_p[j]
-                                       - stokes_phi_p[i] * stokes_div_phi_u[j])
-                                      * fe_values.JxW(q);
-           }
-       }
-      else
-       {
-         const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
-         Assert (dofs_per_cell == elasticity_dofs_per_cell,
-                 ExcInternalError());
+           for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+             {
+               for (unsigned int k=0; k<dofs_per_cell; ++k)
+                 {
+                   elasticity_phi_grad[k] = fe_values[displacements].gradient (k, q);
+                   elasticity_phi_div[k]  = fe_values[displacements].divergence (k, q);
+                 }
 
-         for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
-           {
-             for (unsigned int k=0; k<dofs_per_cell; ++k)
-               {
-                 elasticity_phi_grad[k] = fe_values[displacements].gradient (k, q);
-                 elasticity_phi_div[k]  = fe_values[displacements].divergence (k, q);
-               }
-
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   {
+                     local_matrix(i,j)
+                       +=  (lambda *
+                            elasticity_phi_div[i] * elasticity_phi_div[j]
+                            +
+                            mu *
+                            scalar_product(elasticity_phi_grad[i], elasticity_phi_grad[j])
+                            +
+                            mu *
+                            scalar_product(elasticity_phi_grad[i], transpose(elasticity_phi_grad[j]))
+                       )
+                       *
+                       fe_values.JxW(q);
+                   }
+             }
+         }
+
+                                        // Once we have the contributions from
+                                        // cell integrals, we copy them into
+                                        // the global matrix (taking care of
+                                        // constraints right away, through the
+                                        // ConstraintMatrix::distribute_local_to_global
+                                        // function). Note that we have not
+                                        // written anything into the
+                                        // <code>local_rhs</code> variable,
+                                        // though we still need to pass it
+                                        // along since the elimination of
+                                        // nonzero boundary values requires the
+                                        // modification of local and
+                                        // consequently also global right hand
+                                        // side values:
+       local_dof_indices.resize (cell->get_fe().dofs_per_cell);
+       cell->get_dof_indices (local_dof_indices);
+       constraints.distribute_local_to_global (local_matrix, local_rhs,
+                                               local_dof_indices,
+                                               system_matrix, system_rhs);
+
+                                        // The more interesting part of this
+                                        // function is where we see about face
+                                        // terms along the interface between
+                                        // the two subdomains. To this end, we
+                                        // first have to make sure that we only
+                                        // assemble them once even though a
+                                        // loop over all faces of all cells
+                                        // would encounter each part of the
+                                        // interface twice. We arbitrarily make
+                                        // the decision that we will only
+                                        // evaluate interface terms if the
+                                        // current cell is part of the solid
+                                        // subdomain and if, consequently, a
+                                        // face is not at the boundary and the
+                                        // potential neighbor behind it is part
+                                        // of the fluid domain. Let's start
+                                        // with these conditions:
+       if (cell_is_in_solid_domain (cell))
+         for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+           if (cell->at_boundary(f) == false)
+             {
+                                                // At this point we know that
+                                                // the current cell is a
+                                                // candidate for integration
+                                                // and that a neighbor behind
+                                                // face <code>f</code>
+                                                // exists. There are now three
+                                                // possibilities:
+                                                //
+                                                // - The neighbor is at the
+                                                //   same refinement level and
+                                                //   has no children.
+                                                // - The neighbor has children.
+                                                // - The neighbor is coarser.
+                                                //
+                                                // In all three cases, we are
+                                                // only interested in it if it
+                                                // is part of the fluid
+                                                // subdomain. So let us start
+                                                // with the first and simplest
+                                                // case: if the neighbor is at
+                                                // the same level, has no
+                                                // children, and is a fluid
+                                                // cell, then the two cells
+                                                // share a boundary that is
+                                                // part of the interface along
+                                                // which we want to integrate
+                                                // interface terms. All we have
+                                                // to do is initialize two
+                                                // FEFaceValues object with the
+                                                // current face and the face of
+                                                // the neighboring cell (note
+                                                // how we find out which face
+                                                // of the neighboring cell
+                                                // borders on the current cell)
+                                                // and pass things off to the
+                                                // function that evaluates the
+                                                // interface terms (the third
+                                                // through fifth arguments to
+                                                // this function provide it
+                                                // with scratch arrays). The
+                                                // result is then again copied
+                                                // into the global matrix,
+                                                // using a function that knows
+                                                // that the DoF indices of rows
+                                                // and columns of the local
+                                                // matrix result from different
+                                                // cells:
+               if ((cell->neighbor(f)->level() == cell->level())
+                   &&
+                   (cell->neighbor(f)->has_children() == false)
+                   &&
+                   cell_is_in_fluid_domain (cell->neighbor(f)))
                  {
-                   local_matrix(i,j)
-                     +=  (lambda *
-                          elasticity_phi_div[i] * elasticity_phi_div[j]
-                          +
-                          mu *
-                          scalar_product(elasticity_phi_grad[i], elasticity_phi_grad[j])
-                          +
-                          mu *
-                          scalar_product(elasticity_phi_grad[i], transpose(elasticity_phi_grad[j]))
-                     )
-                     *
-                     fe_values.JxW(q);
+                   elasticity_fe_face_values.reinit (cell, f);
+                   stokes_fe_face_values.reinit (cell->neighbor(f),
+                                                 cell->neighbor_of_neighbor(f));
+
+                   assemble_interface_term (elasticity_fe_face_values, stokes_fe_face_values,
+                                            elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
+                                            local_interface_matrix);
+
+                   cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
+                   constraints.distribute_local_to_global(local_interface_matrix,
+                                                          local_dof_indices,
+                                                          neighbor_dof_indices,
+                                                          system_matrix);
                  }
-           }
-       }
 
-                                      // Once we have the contributions from
-                                      // cell integrals, we copy them into
-                                      // the global matrix (taking care of
-                                      // constraints right away, through the
-                                      // ConstraintMatrix::distribute_local_to_global
-                                      // function). Note that we have not
-                                      // written anything into the
-                                      // <code>local_rhs</code> variable,
-                                      // though we still need to pass it
-                                      // along since the elimination of
-                                      // nonzero boundary values requires the
-                                      // modification of local and
-                                      // consequently also global right hand
-                                      // side values:
-      local_dof_indices.resize (cell->get_fe().dofs_per_cell);
-      cell->get_dof_indices (local_dof_indices);
-      constraints.distribute_local_to_global (local_matrix, local_rhs,
-                                             local_dof_indices,
-                                             system_matrix, system_rhs);
-
-                                      // The more interesting part of this
-                                      // function is where we see about face
-                                      // terms along the interface between
-                                      // the two subdomains. To this end, we
-                                      // first have to make sure that we only
-                                      // assemble them once even though a
-                                      // loop over all faces of all cells
-                                      // would encounter each part of the
-                                      // interface twice. We arbitrarily make
-                                      // the decision that we will only
-                                      // evaluate interface terms if the
-                                      // current cell is part of the solid
-                                      // subdomain and if, consequently, a
-                                      // face is not at the boundary and the
-                                      // potential neighbor behind it is part
-                                      // of the fluid domain. Let's start
-                                      // with these conditions:
-      if (cell_is_in_solid_domain (cell))
-       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-         if (cell->at_boundary(f) == false)
-           {
-                                              // At this point we know that
-                                              // the current cell is a
-                                              // candidate for integration
-                                              // and that a neighbor behind
-                                              // face <code>f</code>
-                                              // exists. There are now three
-                                              // possibilities:
-                                              //
-                                              // - The neighbor is at the
-                                              //   same refinement level and
-                                              //   has no children.
-                                              // - The neighbor has children.
-                                              // - The neighbor is coarser.
-                                              //
-                                              // In all three cases, we are
-                                              // only interested in it if it
-                                              // is part of the fluid
-                                              // subdomain. So let us start
-                                              // with the first and simplest
-                                              // case: if the neighbor is at
-                                              // the same level, has no
-                                              // children, and is a fluid
-                                              // cell, then the two cells
-                                              // share a boundary that is
-                                              // part of the interface along
-                                              // which we want to integrate
-                                              // interface terms. All we have
-                                              // to do is initialize two
-                                              // FEFaceValues object with the
-                                              // current face and the face of
-                                              // the neighboring cell (note
-                                              // how we find out which face
-                                              // of the neighboring cell
-                                              // borders on the current cell)
-                                              // and pass things off to the
-                                              // function that evaluates the
-                                              // interface terms (the third
-                                              // through fifth arguments to
-                                              // this function provide it
-                                              // with scratch arrays). The
-                                              // result is then again copied
-                                              // into the global matrix,
-                                              // using a function that knows
-                                              // that the DoF indices of rows
-                                              // and columns of the local
-                                              // matrix result from different
-                                              // cells:
-             if ((cell->neighbor(f)->level() == cell->level())
-                 &&
-                 (cell->neighbor(f)->has_children() == false)
-                 &&
-                 cell_is_in_fluid_domain (cell->neighbor(f)))
-               {
-                 elasticity_fe_face_values.reinit (cell, f);
-                 stokes_fe_face_values.reinit (cell->neighbor(f),
-                                               cell->neighbor_of_neighbor(f));
-
-                 assemble_interface_term (elasticity_fe_face_values, stokes_fe_face_values,
-                                          elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
-                                          local_interface_matrix);
-
-                 cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
-                 constraints.distribute_local_to_global(local_interface_matrix,
-                                                        local_dof_indices,
-                                                        neighbor_dof_indices,
-                                                        system_matrix);
-               }
-
-                                              // The second case is if the
-                                              // neighbor has further
-                                              // children. In that case, we
-                                              // have to loop over all the
-                                              // children of the neighbor to
-                                              // see if they are part of the
-                                              // fluid subdomain. If they
-                                              // are, then we integrate over
-                                              // the common interface, which
-                                              // is a face for the neighbor
-                                              // and a subface of the current
-                                              // cell, requiring us to use an
-                                              // FEFaceValues for the
-                                              // neighbor and an
-                                              // FESubfaceValues for the
-                                              // current cell:
-             else if ((cell->neighbor(f)->level() == cell->level())
-                      &&
-                      (cell->neighbor(f)->has_children() == true))
-               {
-                 for (unsigned int subface=0;
-                      subface<cell->face(f)->n_children();
-                      ++subface)
-                   if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
-                                                (f, subface)))
-                     {
-                       elasticity_fe_subface_values.reinit (cell,
-                                                            f,
-                                                            subface);
-                       stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface),
-                                                     cell->neighbor_of_neighbor(f));
-
-                       assemble_interface_term (elasticity_fe_subface_values,
-                                                stokes_fe_face_values,
-                                                elasticity_phi,
-                                                stokes_phi_grads_u, stokes_phi_p,
-                                                local_interface_matrix);
-
-                       cell->neighbor_child_on_subface (f, subface)
-                         ->get_dof_indices (neighbor_dof_indices);
-                       constraints.distribute_local_to_global(local_interface_matrix,
-                                                              local_dof_indices,
-                                                              neighbor_dof_indices,
-                                                              system_matrix);
-                     }
-               }
-
-                                              // The last option is that the
-                                              // neighbor is coarser. In that
-                                              // case we have to use an
-                                              // FESubfaceValues object for
-                                              // the neighbor and a
-                                              // FEFaceValues for the current
-                                              // cell; the rest is the same
-                                              // as before:
-             else if (cell->neighbor_is_coarser(f)
-                      &&
-                      cell_is_in_fluid_domain(cell->neighbor(f)))
-               {
-                 elasticity_fe_face_values.reinit (cell, f);
-                 stokes_fe_subface_values.reinit (cell->neighbor(f),
-                                                  cell->neighbor_of_coarser_neighbor(f).first,
-                                                  cell->neighbor_of_coarser_neighbor(f).second);
-
-                 assemble_interface_term (elasticity_fe_face_values,
-                                          stokes_fe_subface_values,
-                                          elasticity_phi,
-                                          stokes_phi_grads_u, stokes_phi_p,
-                                          local_interface_matrix);
-
-                 cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
-                 constraints.distribute_local_to_global(local_interface_matrix,
-                                                        local_dof_indices,
-                                                        neighbor_dof_indices,
-                                                        system_matrix);
-
-               }
-           }
-    }
-}
+                                                // The second case is if the
+                                                // neighbor has further
+                                                // children. In that case, we
+                                                // have to loop over all the
+                                                // children of the neighbor to
+                                                // see if they are part of the
+                                                // fluid subdomain. If they
+                                                // are, then we integrate over
+                                                // the common interface, which
+                                                // is a face for the neighbor
+                                                // and a subface of the current
+                                                // cell, requiring us to use an
+                                                // FEFaceValues for the
+                                                // neighbor and an
+                                                // FESubfaceValues for the
+                                                // current cell:
+               else if ((cell->neighbor(f)->level() == cell->level())
+                        &&
+                        (cell->neighbor(f)->has_children() == true))
+                 {
+                   for (unsigned int subface=0;
+                        subface<cell->face(f)->n_children();
+                        ++subface)
+                     if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+                                                  (f, subface)))
+                       {
+                         elasticity_fe_subface_values.reinit (cell,
+                                                              f,
+                                                              subface);
+                         stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface),
+                                                       cell->neighbor_of_neighbor(f));
+
+                         assemble_interface_term (elasticity_fe_subface_values,
+                                                  stokes_fe_face_values,
+                                                  elasticity_phi,
+                                                  stokes_phi_grads_u, stokes_phi_p,
+                                                  local_interface_matrix);
+
+                         cell->neighbor_child_on_subface (f, subface)
+                           ->get_dof_indices (neighbor_dof_indices);
+                         constraints.distribute_local_to_global(local_interface_matrix,
+                                                                local_dof_indices,
+                                                                neighbor_dof_indices,
+                                                                system_matrix);
+                       }
+                 }
 
+                                                // The last option is that the
+                                                // neighbor is coarser. In that
+                                                // case we have to use an
+                                                // FESubfaceValues object for
+                                                // the neighbor and a
+                                                // FEFaceValues for the current
+                                                // cell; the rest is the same
+                                                // as before:
+               else if (cell->neighbor_is_coarser(f)
+                        &&
+                        cell_is_in_fluid_domain(cell->neighbor(f)))
+                 {
+                   elasticity_fe_face_values.reinit (cell, f);
+                   stokes_fe_subface_values.reinit (cell->neighbor(f),
+                                                    cell->neighbor_of_coarser_neighbor(f).first,
+                                                    cell->neighbor_of_coarser_neighbor(f).second);
+
+                   assemble_interface_term (elasticity_fe_face_values,
+                                            stokes_fe_subface_values,
+                                            elasticity_phi,
+                                            stokes_phi_grads_u, stokes_phi_p,
+                                            local_interface_matrix);
+
+                   cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
+                   constraints.distribute_local_to_global(local_interface_matrix,
+                                                          local_dof_indices,
+                                                          neighbor_dof_indices,
+                                                          system_matrix);
 
+                 }
+             }
+      }
+  }
 
-                                // In the function that assembles the global
-                                // system, we passed computing interface
-                                // terms to a separate function we discuss
-                                // here. The key is that even though we can't
-                                // predict the combination of FEFaceValues
-                                // and FESubfaceValues objects, they are both
-                                // derived from the FEFaceValuesBase class
-                                // and consequently we don't have to care:
-                                // the function is simply called with two
-                                // such objects denoting the values of the
-                                // shape functions on the quadrature points
-                                // of the two sides of the face. We then do
-                                // what we always do: we fill the scratch
-                                // arrays with the values of shape functions
-                                // and their derivatives, and then loop over
-                                // all entries of the matrix to compute the
-                                // local integrals. The details of the
-                                // bilinear form we evaluate here are given
-                                // in the introduction.
-template <int dim>
-void
-FluidStructureProblem<dim>::
-assemble_interface_term (const FEFaceValuesBase<dim>          &elasticity_fe_face_values,
-                        const FEFaceValuesBase<dim>          &stokes_fe_face_values,
-                        std::vector<Tensor<1,dim> >          &elasticity_phi,
-                        std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
-                        std::vector<double>                  &stokes_phi_p,
-                        FullMatrix<double>                   &local_interface_matrix) const
-{
-  Assert (stokes_fe_face_values.n_quadrature_points ==
-          elasticity_fe_face_values.n_quadrature_points,
-         ExcInternalError());
-  const unsigned int n_face_quadrature_points
-    = elasticity_fe_face_values.n_quadrature_points;
-
-  const FEValuesExtractors::Vector velocities (0);
-  const FEValuesExtractors::Scalar pressure (dim);
-  const FEValuesExtractors::Vector displacements (dim+1);
-
-  local_interface_matrix = 0;
-  for (unsigned int q=0; q<n_face_quadrature_points; ++q)
-    {
-      const Tensor<1,dim> normal_vector = stokes_fe_face_values.normal_vector(q);
-
-      for (unsigned int k=0; k<stokes_fe_face_values.dofs_per_cell; ++k)
-       stokes_phi_grads_u[k] = stokes_fe_face_values[velocities].symmetric_gradient (k, q);
-      for (unsigned int k=0; k<elasticity_fe_face_values.dofs_per_cell; ++k)
-       elasticity_phi[k] = elasticity_fe_face_values[displacements].value (k,q);
-
-      for (unsigned int i=0; i<elasticity_fe_face_values.dofs_per_cell; ++i)
-       for (unsigned int j=0; j<stokes_fe_face_values.dofs_per_cell; ++j)
-         local_interface_matrix(i,j) += -((2 * viscosity *
-                                           (stokes_phi_grads_u[j] *
-                                            normal_vector)
-                                           +
-                                           stokes_phi_p[j] *
-                                           normal_vector) *
-                                          elasticity_phi[i] *
-                                          stokes_fe_face_values.JxW(q));
-    }
-}
 
 
-                                 // @sect4{<code>FluidStructureProblem::solve</code>}
+                                  // In the function that assembles the global
+                                  // system, we passed computing interface
+                                  // terms to a separate function we discuss
+                                  // here. The key is that even though we can't
+                                  // predict the combination of FEFaceValues
+                                  // and FESubfaceValues objects, they are both
+                                  // derived from the FEFaceValuesBase class
+                                  // and consequently we don't have to care:
+                                  // the function is simply called with two
+                                  // such objects denoting the values of the
+                                  // shape functions on the quadrature points
+                                  // of the two sides of the face. We then do
+                                  // what we always do: we fill the scratch
+                                  // arrays with the values of shape functions
+                                  // and their derivatives, and then loop over
+                                  // all entries of the matrix to compute the
+                                  // local integrals. The details of the
+                                  // bilinear form we evaluate here are given
+                                  // in the introduction.
+  template <int dim>
+  void
+  FluidStructureProblem<dim>::
+  assemble_interface_term (const FEFaceValuesBase<dim>          &elasticity_fe_face_values,
+                          const FEFaceValuesBase<dim>          &stokes_fe_face_values,
+                          std::vector<Tensor<1,dim> >          &elasticity_phi,
+                          std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
+                          std::vector<double>                  &stokes_phi_p,
+                          FullMatrix<double>                   &local_interface_matrix) const
+  {
+    Assert (stokes_fe_face_values.n_quadrature_points ==
+           elasticity_fe_face_values.n_quadrature_points,
+           ExcInternalError());
+    const unsigned int n_face_quadrature_points
+      = elasticity_fe_face_values.n_quadrature_points;
+
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
+    const FEValuesExtractors::Vector displacements (dim+1);
+
+    local_interface_matrix = 0;
+    for (unsigned int q=0; q<n_face_quadrature_points; ++q)
+      {
+       const Tensor<1,dim> normal_vector = stokes_fe_face_values.normal_vector(q);
+
+       for (unsigned int k=0; k<stokes_fe_face_values.dofs_per_cell; ++k)
+         stokes_phi_grads_u[k] = stokes_fe_face_values[velocities].symmetric_gradient (k, q);
+       for (unsigned int k=0; k<elasticity_fe_face_values.dofs_per_cell; ++k)
+         elasticity_phi[k] = elasticity_fe_face_values[displacements].value (k,q);
+
+       for (unsigned int i=0; i<elasticity_fe_face_values.dofs_per_cell; ++i)
+         for (unsigned int j=0; j<stokes_fe_face_values.dofs_per_cell; ++j)
+           local_interface_matrix(i,j) += -((2 * viscosity *
+                                             (stokes_phi_grads_u[j] *
+                                              normal_vector)
+                                             +
+                                             stokes_phi_p[j] *
+                                             normal_vector) *
+                                            elasticity_phi[i] *
+                                            stokes_fe_face_values.JxW(q));
+      }
+  }
 
-                                // As discussed in the introduction, we use a
-                                // rather trivial solver here: we just pass
-                                // the linear system off to the
-                                // SparseDirectUMFPACK direct solver (see,
-                                // for example, step-29). The only thing we
-                                // have to do after solving is ensure that
-                                // hanging node and boundary value
-                                // constraints are correct.
-template <int dim>
-void
-FluidStructureProblem<dim>::solve ()
-{
-  SparseDirectUMFPACK direct_solver;
-  direct_solver.initialize (system_matrix);
-  direct_solver.vmult (solution, system_rhs);
 
-  constraints.distribute (solution);
-}
+                                  // @sect4{<code>FluidStructureProblem::solve</code>}
+
+                                  // As discussed in the introduction, we use a
+                                  // rather trivial solver here: we just pass
+                                  // the linear system off to the
+                                  // SparseDirectUMFPACK direct solver (see,
+                                  // for example, step-29). The only thing we
+                                  // have to do after solving is ensure that
+                                  // hanging node and boundary value
+                                  // constraints are correct.
+  template <int dim>
+  void
+  FluidStructureProblem<dim>::solve ()
+  {
+    SparseDirectUMFPACK direct_solver;
+    direct_solver.initialize (system_matrix);
+    direct_solver.vmult (solution, system_rhs);
 
+    constraints.distribute (solution);
+  }
 
 
-                                 // @sect4{<code>FluidStructureProblem::output_results</code>}
 
-                                // Generating graphical output is rather
-                                // trivial here: all we have to do is
-                                // identify which components of the solution
-                                // vector belong to scalars and/or vectors
-                                // (see, for example, step-22 for a previous
-                                // example), and then pass it all on to the
-                                // DataOut class (with the second template
-                                // argument equal to hp::DoFHandler instead
-                                // of the usual default DoFHandler):
-template <int dim>
-void
-FluidStructureProblem<dim>::
-output_results (const unsigned int refinement_cycle)  const
-{
-  std::vector<std::string> solution_names (dim, "velocity");
-  solution_names.push_back ("pressure");
-  for (unsigned int d=0; d<dim; ++d)
-    solution_names.push_back ("displacement");
+                                  // @sect4{<code>FluidStructureProblem::output_results</code>}
 
-  std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    data_component_interpretation
-    (dim, DataComponentInterpretation::component_is_part_of_vector);
-  data_component_interpretation
-    .push_back (DataComponentInterpretation::component_is_scalar);
-  for (unsigned int d=0; d<dim; ++d)
+                                  // Generating graphical output is rather
+                                  // trivial here: all we have to do is
+                                  // identify which components of the solution
+                                  // vector belong to scalars and/or vectors
+                                  // (see, for example, step-22 for a previous
+                                  // example), and then pass it all on to the
+                                  // DataOut class (with the second template
+                                  // argument equal to hp::DoFHandler instead
+                                  // of the usual default DoFHandler):
+  template <int dim>
+  void
+  FluidStructureProblem<dim>::
+  output_results (const unsigned int refinement_cycle)  const
+  {
+    std::vector<std::string> solution_names (dim, "velocity");
+    solution_names.push_back ("pressure");
+    for (unsigned int d=0; d<dim; ++d)
+      solution_names.push_back ("displacement");
+
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      data_component_interpretation
+      (dim, DataComponentInterpretation::component_is_part_of_vector);
     data_component_interpretation
-      .push_back (DataComponentInterpretation::component_is_part_of_vector);
+      .push_back (DataComponentInterpretation::component_is_scalar);
+    for (unsigned int d=0; d<dim; ++d)
+      data_component_interpretation
+       .push_back (DataComponentInterpretation::component_is_part_of_vector);
 
-  DataOut<dim,hp::DoFHandler<dim> > data_out;
-  data_out.attach_dof_handler (dof_handler);
+    DataOut<dim,hp::DoFHandler<dim> > data_out;
+    data_out.attach_dof_handler (dof_handler);
 
-  data_out.add_data_vector (solution, solution_names,
-                           DataOut<dim,hp::DoFHandler<dim> >::type_dof_data,
-                           data_component_interpretation);
-  data_out.build_patches ();
+    data_out.add_data_vector (solution, solution_names,
+                             DataOut<dim,hp::DoFHandler<dim> >::type_dof_data,
+                             data_component_interpretation);
+    data_out.build_patches ();
 
-  std::ostringstream filename;
-  filename << "solution-"
-           << Utilities::int_to_string (refinement_cycle, 2)
-           << ".vtk";
+    std::ostringstream filename;
+    filename << "solution-"
+            << Utilities::int_to_string (refinement_cycle, 2)
+            << ".vtk";
 
-  std::ofstream output (filename.str().c_str());
-  data_out.write_vtk (output);
-}
+    std::ofstream output (filename.str().c_str());
+    data_out.write_vtk (output);
+  }
 
 
-                                 // @sect4{<code>FluidStructureProblem::refine_mesh</code>}
-
-                                // The next step is to refine the mesh. As
-                                // was discussed in the introduction, this is
-                                // a bit tricky primarily because the fluid
-                                // and the solid subdomains use variables
-                                // that have different physical dimensions
-                                // and for which the absolute magnitude of
-                                // error estimates is consequently not
-                                // directly comparable. We will therefore
-                                // have to scale them. At the top of the
-                                // function, we therefore first compute error
-                                // estimates for the different variables
-                                // separately (using the velocities but not
-                                // the pressure for the fluid domain, and the
-                                // displacements in the solid domain):
-template <int dim>
-void
-FluidStructureProblem<dim>::refine_mesh ()
-{
-  Vector<float>
-    stokes_estimated_error_per_cell (triangulation.n_active_cells());
-  Vector<float>
-    elasticity_estimated_error_per_cell (triangulation.n_active_cells());
-
-  const QGauss<dim-1> stokes_face_quadrature(stokes_degree+2);
-  const QGauss<dim-1> elasticity_face_quadrature(elasticity_degree+2);
-
-  hp::QCollection<dim-1> face_q_collection;
-  face_q_collection.push_back (stokes_face_quadrature);
-  face_q_collection.push_back (elasticity_face_quadrature);
-
-  std::vector<bool> stokes_component_mask (dim+1+dim, false);
-  for (unsigned int d=0; d<dim; ++d)
-    stokes_component_mask[d] = true;
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                      face_q_collection,
-                                      typename FunctionMap<dim>::type(),
-                                      solution,
-                                      stokes_estimated_error_per_cell,
-                                      stokes_component_mask);
-
-  std::vector<bool> elasticity_component_mask (dim+1+dim, false);
-  for (unsigned int d=0; d<dim; ++d)
-    elasticity_component_mask[dim+1+d] = true;
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                      face_q_collection,
-                                      typename FunctionMap<dim>::type(),
-                                      solution,
-                                      elasticity_estimated_error_per_cell,
-                                      elasticity_component_mask);
-
-                                  // We then normalize error estimates by
-                                  // dividing by their norm and scale the
-                                  // fluid error indicators by a factor of 4
-                                  // as discussed in the introduction. The
-                                  // results are then added together into a
-                                  // vector that contains error indicators
-                                  // for all cells:
-  stokes_estimated_error_per_cell
-    *= 4. / stokes_estimated_error_per_cell.l2_norm();
-  elasticity_estimated_error_per_cell
-    *= 1. / elasticity_estimated_error_per_cell.l2_norm();
-
-  Vector<float>
-    estimated_error_per_cell (triangulation.n_active_cells());
-
-  estimated_error_per_cell += stokes_estimated_error_per_cell;
-  estimated_error_per_cell += elasticity_estimated_error_per_cell;
-
-                                  // The second to last part of the function,
-                                  // before actually refining the mesh,
-                                  // involves a heuristic that we have
-                                  // already mentioned in the introduction:
-                                  // because the solution is discontinuous,
-                                  // the KellyErrorEstimator class gets all
-                                  // confused about cells that sit at the
-                                  // boundary between subdomains: it believes
-                                  // that the error is large there because
-                                  // the jump in the gradient is large, even
-                                  // though this is entirely expected and a
-                                  // feature that is in fact present in the
-                                  // exact solution as well and therefore not
-                                  // indicative of any numerical error.
-                                  //
-                                  // Consequently, we set the error
-                                  // indicators to zero for all cells at the
-                                  // interface; the conditions determining
-                                  // which cells this affects are slightly
-                                  // awkward because we have to account for
-                                  // the possibility of adaptively refined
-                                  // meshes, meaning that the neighboring
-                                  // cell can be coarser than the current
-                                  // one, or could in fact be refined some
-                                  // more. The structure of these nested
-                                  // conditions is much the same as we
-                                  // encountered when assembling interface
-                                  // terms in <code>assemble_system</code>.
+                                  // @sect4{<code>FluidStructureProblem::refine_mesh</code>}
+
+                                  // The next step is to refine the mesh. As
+                                  // was discussed in the introduction, this is
+                                  // a bit tricky primarily because the fluid
+                                  // and the solid subdomains use variables
+                                  // that have different physical dimensions
+                                  // and for which the absolute magnitude of
+                                  // error estimates is consequently not
+                                  // directly comparable. We will therefore
+                                  // have to scale them. At the top of the
+                                  // function, we therefore first compute error
+                                  // estimates for the different variables
+                                  // separately (using the velocities but not
+                                  // the pressure for the fluid domain, and the
+                                  // displacements in the solid domain):
+  template <int dim>
+  void
+  FluidStructureProblem<dim>::refine_mesh ()
   {
-    unsigned int cell_index = 0;
-    for (typename hp::DoFHandler<dim>::active_cell_iterator
-          cell = dof_handler.begin_active();
-        cell != dof_handler.end(); ++cell, ++cell_index)
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-       if (cell_is_in_solid_domain (cell))
-         {
-           if ((cell->at_boundary(f) == false)
-               &&
-               (((cell->neighbor(f)->level() == cell->level())
-                 &&
-                 (cell->neighbor(f)->has_children() == false)
-                 &&
-                 cell_is_in_fluid_domain (cell->neighbor(f)))
-                ||
-                ((cell->neighbor(f)->level() == cell->level())
-                 &&
-                 (cell->neighbor(f)->has_children() == true)
-                 &&
-                 (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
-                                           (f, 0))))
-                ||
-                (cell->neighbor_is_coarser(f)
-                 &&
-                 cell_is_in_fluid_domain(cell->neighbor(f)))
-               ))
-             estimated_error_per_cell(cell_index) = 0;
-         }
-       else
-         {
-           if ((cell->at_boundary(f) == false)
-               &&
-               (((cell->neighbor(f)->level() == cell->level())
-                 &&
-                 (cell->neighbor(f)->has_children() == false)
-                 &&
-                 cell_is_in_solid_domain (cell->neighbor(f)))
-                ||
-                ((cell->neighbor(f)->level() == cell->level())
-                 &&
-                 (cell->neighbor(f)->has_children() == true)
+    Vector<float>
+      stokes_estimated_error_per_cell (triangulation.n_active_cells());
+    Vector<float>
+      elasticity_estimated_error_per_cell (triangulation.n_active_cells());
+
+    const QGauss<dim-1> stokes_face_quadrature(stokes_degree+2);
+    const QGauss<dim-1> elasticity_face_quadrature(elasticity_degree+2);
+
+    hp::QCollection<dim-1> face_q_collection;
+    face_q_collection.push_back (stokes_face_quadrature);
+    face_q_collection.push_back (elasticity_face_quadrature);
+
+    std::vector<bool> stokes_component_mask (dim+1+dim, false);
+    for (unsigned int d=0; d<dim; ++d)
+      stokes_component_mask[d] = true;
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                       face_q_collection,
+                                       typename FunctionMap<dim>::type(),
+                                       solution,
+                                       stokes_estimated_error_per_cell,
+                                       stokes_component_mask);
+
+    std::vector<bool> elasticity_component_mask (dim+1+dim, false);
+    for (unsigned int d=0; d<dim; ++d)
+      elasticity_component_mask[dim+1+d] = true;
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                       face_q_collection,
+                                       typename FunctionMap<dim>::type(),
+                                       solution,
+                                       elasticity_estimated_error_per_cell,
+                                       elasticity_component_mask);
+
+                                    // We then normalize error estimates by
+                                    // dividing by their norm and scale the
+                                    // fluid error indicators by a factor of 4
+                                    // as discussed in the introduction. The
+                                    // results are then added together into a
+                                    // vector that contains error indicators
+                                    // for all cells:
+    stokes_estimated_error_per_cell
+      *= 4. / stokes_estimated_error_per_cell.l2_norm();
+    elasticity_estimated_error_per_cell
+      *= 1. / elasticity_estimated_error_per_cell.l2_norm();
+
+    Vector<float>
+      estimated_error_per_cell (triangulation.n_active_cells());
+
+    estimated_error_per_cell += stokes_estimated_error_per_cell;
+    estimated_error_per_cell += elasticity_estimated_error_per_cell;
+
+                                    // The second to last part of the function,
+                                    // before actually refining the mesh,
+                                    // involves a heuristic that we have
+                                    // already mentioned in the introduction:
+                                    // because the solution is discontinuous,
+                                    // the KellyErrorEstimator class gets all
+                                    // confused about cells that sit at the
+                                    // boundary between subdomains: it believes
+                                    // that the error is large there because
+                                    // the jump in the gradient is large, even
+                                    // though this is entirely expected and a
+                                    // feature that is in fact present in the
+                                    // exact solution as well and therefore not
+                                    // indicative of any numerical error.
+                                    //
+                                    // Consequently, we set the error
+                                    // indicators to zero for all cells at the
+                                    // interface; the conditions determining
+                                    // which cells this affects are slightly
+                                    // awkward because we have to account for
+                                    // the possibility of adaptively refined
+                                    // meshes, meaning that the neighboring
+                                    // cell can be coarser than the current
+                                    // one, or could in fact be refined some
+                                    // more. The structure of these nested
+                                    // conditions is much the same as we
+                                    // encountered when assembling interface
+                                    // terms in <code>assemble_system</code>.
+    {
+      unsigned int cell_index = 0;
+      for (typename hp::DoFHandler<dim>::active_cell_iterator
+            cell = dof_handler.begin_active();
+          cell != dof_handler.end(); ++cell, ++cell_index)
+       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+         if (cell_is_in_solid_domain (cell))
+           {
+             if ((cell->at_boundary(f) == false)
                  &&
-                 (cell_is_in_solid_domain (cell->neighbor_child_on_subface
-                                           (f, 0))))
-                ||
-                (cell->neighbor_is_coarser(f)
+                 (((cell->neighbor(f)->level() == cell->level())
+                   &&
+                   (cell->neighbor(f)->has_children() == false)
+                   &&
+                   cell_is_in_fluid_domain (cell->neighbor(f)))
+                  ||
+                  ((cell->neighbor(f)->level() == cell->level())
+                   &&
+                   (cell->neighbor(f)->has_children() == true)
+                   &&
+                   (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+                                             (f, 0))))
+                  ||
+                  (cell->neighbor_is_coarser(f)
+                   &&
+                   cell_is_in_fluid_domain(cell->neighbor(f)))
+                 ))
+               estimated_error_per_cell(cell_index) = 0;
+           }
+         else
+           {
+             if ((cell->at_boundary(f) == false)
                  &&
-                 cell_is_in_solid_domain(cell->neighbor(f)))
-               ))
-             estimated_error_per_cell(cell_index) = 0;
-         }
-  }
+                 (((cell->neighbor(f)->level() == cell->level())
+                   &&
+                   (cell->neighbor(f)->has_children() == false)
+                   &&
+                   cell_is_in_solid_domain (cell->neighbor(f)))
+                  ||
+                  ((cell->neighbor(f)->level() == cell->level())
+                   &&
+                   (cell->neighbor(f)->has_children() == true)
+                   &&
+                   (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+                                             (f, 0))))
+                  ||
+                  (cell->neighbor_is_coarser(f)
+                   &&
+                   cell_is_in_solid_domain(cell->neighbor(f)))
+                 ))
+               estimated_error_per_cell(cell_index) = 0;
+           }
+    }
 
-  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                   estimated_error_per_cell,
-                                                   0.3, 0.0);
-  triangulation.execute_coarsening_and_refinement ();
-}
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.3, 0.0);
+    triangulation.execute_coarsening_and_refinement ();
+  }
 
 
 
-                                 // @sect4{<code>FluidStructureProblem::run</code>}
+                                  // @sect4{<code>FluidStructureProblem::run</code>}
 
-                                // This is, as usual, the function that
-                                // controls the overall flow of operation. If
-                                // you've read through tutorial programs
-                                // step-1 through step-6, for example, then
-                                // you are already quite familiar with the
-                                // following structure:
-template <int dim>
-void FluidStructureProblem<dim>::run ()
-{
-  make_grid ();
+                                  // This is, as usual, the function that
+                                  // controls the overall flow of operation. If
+                                  // you've read through tutorial programs
+                                  // step-1 through step-6, for example, then
+                                  // you are already quite familiar with the
+                                  // following structure:
+  template <int dim>
+  void FluidStructureProblem<dim>::run ()
+  {
+    make_grid ();
 
-  for (unsigned int refinement_cycle = 0; refinement_cycle<10-2*dim;
-       ++refinement_cycle)
-    {
-      std::cout << "Refinement cycle " << refinement_cycle << std::endl;
+    for (unsigned int refinement_cycle = 0; refinement_cycle<10-2*dim;
+        ++refinement_cycle)
+      {
+       std::cout << "Refinement cycle " << refinement_cycle << std::endl;
 
-      if (refinement_cycle > 0)
-        refine_mesh ();
+       if (refinement_cycle > 0)
+         refine_mesh ();
 
-      setup_dofs ();
+       setup_dofs ();
 
-      std::cout << "   Assembling..." << std::endl;
-      assemble_system ();
+       std::cout << "   Assembling..." << std::endl;
+       assemble_system ();
 
-      std::cout << "   Solving..." << std::endl;
-      solve ();
+       std::cout << "   Solving..." << std::endl;
+       solve ();
 
-      std::cout << "   Writing output..." << std::endl;
-      output_results (refinement_cycle);
+       std::cout << "   Writing output..." << std::endl;
+       output_results (refinement_cycle);
 
-      std::cout << std::endl;
-    }
+       std::cout << std::endl;
+      }
+  }
 }
 
 
@@ -1259,6 +1263,9 @@ int main ()
 {
   try
     {
+      using namespace dealii;
+      using namespace Step46;
+
       deallog.depth_console (0);
 
       FluidStructureProblem<2> flow_problem(1, 1);
index a442fb28bafa83b46c832c664cba38a327bfff94..0580e15af93442517170a1b3bff7fe4c8ad7a6b9 100644 (file)
 
 #include <deal.II/numerics/error_estimator.h>
 
-using namespace dealii;
-
-
-
-double sign (double d)
-{
-  if (d > 0)
-    return 1;
-  else if (d < 0)
-    return -1;
-  else
-    return 0;
-}
-
-
-template <int dim>
-class LaplaceProblem
+namespace Step47
 {
-  public:
-    LaplaceProblem ();
-    ~LaplaceProblem ();
-
-    void run ();
+  using namespace dealii;
 
-  private:
-    bool interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const;
-    std::pair<unsigned int, Quadrature<dim> > compute_quadrature(const Quadrature<dim> &plain_quadrature, const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const std::vector<double> &level_set_values);
-    void append_quadrature(const Quadrature<dim> &plain_quadrature,
-                          const std::vector<Point<dim> > &v      ,
-                          std::vector<Point<dim> > &xfem_points,
-                          std::vector<double>      &xfem_weights);
 
-    void setup_system ();
-    void assemble_system ();
-    void solve ();
-    void refine_grid ();
-    void output_results (const unsigned int cycle) const;
-    void compute_error () const;
 
-    Triangulation<dim>    triangulation;
+  double sign (double d)
+  {
+    if (d > 0)
+      return 1;
+    else if (d < 0)
+      return -1;
+    else
+      return 0;
+  }
 
-    hp::DoFHandler<dim>   dof_handler;
-    hp::FECollection<dim> fe_collection;
 
-    ConstraintMatrix      constraints;
+  template <int dim>
+  class LaplaceProblem
+  {
+    public:
+      LaplaceProblem ();
+      ~LaplaceProblem ();
 
-    SparsityPattern       sparsity_pattern;
-    SparseMatrix<double>  system_matrix;
+      void run ();
 
-    Vector<double>        solution;
-    Vector<double>        system_rhs;
-};
+    private:
+      bool interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const;
+      std::pair<unsigned int, Quadrature<dim> > compute_quadrature(const Quadrature<dim> &plain_quadrature, const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const std::vector<double> &level_set_values);
+      void append_quadrature(const Quadrature<dim> &plain_quadrature,
+                            const std::vector<Point<dim> > &v      ,
+                            std::vector<Point<dim> > &xfem_points,
+                            std::vector<double>      &xfem_weights);
 
+      void setup_system ();
+      void assemble_system ();
+      void solve ();
+      void refine_grid ();
+      void output_results (const unsigned int cycle) const;
+      void compute_error () const;
 
+      Triangulation<dim>    triangulation;
 
+      hp::DoFHandler<dim>   dof_handler;
+      hp::FECollection<dim> fe_collection;
 
-template <int dim>
-class Coefficient : public Function<dim>
-{
-  public:
-    Coefficient () : Function<dim>() {}
+      ConstraintMatrix      constraints;
 
-    virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
+      SparsityPattern       sparsity_pattern;
+      SparseMatrix<double>  system_matrix;
 
-    virtual void value_list (const std::vector<Point<dim> > &points,
-                            std::vector<double>            &values,
-                            const unsigned int              component = 0) const;
-};
+      Vector<double>        solution;
+      Vector<double>        system_rhs;
+  };
 
 
 
-template <int dim>
-double Coefficient<dim>::value (const Point<dim> &p,
-                               const unsigned int) const
-{
-  if (p.square() < 0.5*0.5)
-    return 20;
-  else
-    return 1;
-}
 
+  template <int dim>
+  class Coefficient : public Function<dim>
+  {
+    public:
+      Coefficient () : Function<dim>() {}
 
+      virtual double value (const Point<dim>   &p,
+                           const unsigned int  component = 0) const;
 
-template <int dim>
-void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
-                                  std::vector<double>            &values,
-                                  const unsigned int              component) const
-{
-  const unsigned int n_points = points.size();
+      virtual void value_list (const std::vector<Point<dim> > &points,
+                              std::vector<double>            &values,
+                              const unsigned int              component = 0) const;
+  };
 
-  Assert (values.size() == n_points,
-         ExcDimensionMismatch (values.size(), n_points));
 
-  Assert (component == 0,
-         ExcIndexRange (component, 0, 1));
 
-  for (unsigned int i=0; i<n_points; ++i)
-    {
-                       if (points[i].square() < 0.5*0.5)
-                               values[i] = 20;
-                       else
-                               values[i] = 1;
-               }
-}
+  template <int dim>
+  double Coefficient<dim>::value (const Point<dim> &p,
+                                 const unsigned int) const
+  {
+    if (p.square() < 0.5*0.5)
+      return 20;
+    else
+      return 1;
+  }
 
 
 
-template <int dim>
-double exact_solution (const Point<dim> &p)
-{
-  const double r = p.norm();
+  template <int dim>
+  void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+                                    std::vector<double>            &values,
+                                    const unsigned int              component) const
+  {
+    const unsigned int n_points = points.size();
 
-  return (r < 0.5
-         ?
-         1./20 * (-1./4*r*r + 61./16)
-         :
-         1./4 * (1-r*r));
-}
+    Assert (values.size() == n_points,
+           ExcDimensionMismatch (values.size(), n_points));
 
+    Assert (component == 0,
+           ExcIndexRange (component, 0, 1));
 
-template <int dim>
-LaplaceProblem<dim>::LaplaceProblem ()
-               :
-               dof_handler (triangulation)
-{
-  fe_collection.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
-                                         FE_Nothing<dim>(), 1));
-  fe_collection.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
-                                         FE_Q<dim>(1), 1));
-}
+    for (unsigned int i=0; i<n_points; ++i)
+      {
+       if (points[i].square() < 0.5*0.5)
+         values[i] = 20;
+       else
+         values[i] = 1;
+      }
+  }
 
 
 
-template <int dim>
-LaplaceProblem<dim>::~LaplaceProblem ()
-{
-  dof_handler.clear ();
-}
+  template <int dim>
+  double exact_solution (const Point<dim> &p)
+  {
+    const double r = p.norm();
 
+    return (r < 0.5
+           ?
+           1./20 * (-1./4*r*r + 61./16)
+           :
+           1./4 * (1-r*r));
+  }
 
 
-template <int dim>
-double
-level_set (const Point<dim> &p)
-{
-  return p.norm() - 0.5;
-}
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem ()
+                 :
+                 dof_handler (triangulation)
+  {
+    fe_collection.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+                                           FE_Nothing<dim>(), 1));
+    fe_collection.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+                                           FE_Q<dim>(1), 1));
+  }
 
 
 
-template <int dim>
-Tensor<1,dim>
-grad_level_set (const Point<dim> &p)
-{
-  return p / p.norm();
-}
+  template <int dim>
+  LaplaceProblem<dim>::~LaplaceProblem ()
+  {
+    dof_handler.clear ();
+  }
 
 
 
-template <int dim>
-bool
-LaplaceProblem<dim>::
-interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const
-{
-  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell-1; ++v)
-    if (level_set(cell->vertex(v)) * level_set(cell->vertex(v+1)) < 0)
-      return true;
-
-                                  // we get here only if all vertices
-                                  // have the same sign, which means
-                                  // that the cell is not intersected
-  return false;
-}
+  template <int dim>
+  double
+  level_set (const Point<dim> &p)
+  {
+    return p.norm() - 0.5;
+  }
 
 
 
-template <int dim>
-void LaplaceProblem<dim>::setup_system ()
-{
-  for (typename hp::DoFHandler<dim>::cell_iterator cell
-        = dof_handler.begin_active();
-       cell != dof_handler.end(); ++cell)
-    if (interface_intersects_cell(cell) == false)
-      cell->set_active_fe_index(0);
-    else
-      cell->set_active_fe_index(1);
+  template <int dim>
+  Tensor<1,dim>
+  grad_level_set (const Point<dim> &p)
+  {
+    return p / p.norm();
+  }
 
-  dof_handler.distribute_dofs (fe_collection);
 
-  solution.reinit (dof_handler.n_dofs());
-  system_rhs.reinit (dof_handler.n_dofs());
 
+  template <int dim>
+  bool
+  LaplaceProblem<dim>::
+  interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const
+  {
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell-1; ++v)
+      if (level_set(cell->vertex(v)) * level_set(cell->vertex(v+1)) < 0)
+       return true;
 
-  constraints.clear ();
-//TODO: fix this, it currently crashes
-  // DoFTools::make_hanging_node_constraints (dof_handler,
-  //                                      constraints);
+                                    // we get here only if all vertices
+                                    // have the same sign, which means
+                                    // that the cell is not intersected
+    return false;
+  }
 
-//TODO: component 1 must satisfy zero boundary conditions
-  constraints.close();
 
 
-  CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
-  DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
+  template <int dim>
+  void LaplaceProblem<dim>::setup_system ()
+  {
+    for (typename hp::DoFHandler<dim>::cell_iterator cell
+          = dof_handler.begin_active();
+        cell != dof_handler.end(); ++cell)
+      if (interface_intersects_cell(cell) == false)
+       cell->set_active_fe_index(0);
+      else
+       cell->set_active_fe_index(1);
 
-  constraints.condense (c_sparsity);
+    dof_handler.distribute_dofs (fe_collection);
 
-  sparsity_pattern.copy_from(c_sparsity);
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
 
-  system_matrix.reinit (sparsity_pattern);
-}
 
+    constraints.clear ();
+//TODO: fix this, it currently crashes
+                                    // DoFTools::make_hanging_node_constraints (dof_handler,
+                                    //                                            constraints);
 
-template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
-  const QGauss<dim>  quadrature_formula(3);
+//TODO: component 1 must satisfy zero boundary conditions
+    constraints.close();
 
 
-  FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
-                                update_values    |  update_gradients |
-                                update_quadrature_points  |  update_JxW_values);
-  FEValues<dim> enriched_fe_values (fe_collection[1], quadrature_formula,
-                                            update_values    |  update_gradients |
-                                            update_quadrature_points  |  update_JxW_values);
+    CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
 
-  const unsigned int   n_q_points    = quadrature_formula.size();
+    constraints.condense (c_sparsity);
 
-  FullMatrix<double>   cell_matrix;
-  Vector<double>       cell_rhs;
+    sparsity_pattern.copy_from(c_sparsity);
 
-  std::vector<unsigned int> local_dof_indices;
+    system_matrix.reinit (sparsity_pattern);
+  }
 
-  const Coefficient<dim> coefficient;
-  std::vector<double>    coefficient_values (n_q_points);
 
-  typename hp::DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
+  template <int dim>
+  void LaplaceProblem<dim>::assemble_system ()
+  {
+    const QGauss<dim>  quadrature_formula(3);
 
-  for (; cell!=endc; ++cell)
-    {
-      const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
-      cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
-      cell_rhs.reinit (dofs_per_cell);
 
-      cell_matrix = 0;
-      cell_rhs = 0;
+    FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
+                                  update_values    |  update_gradients |
+                                  update_quadrature_points  |  update_JxW_values);
+    FEValues<dim> enriched_fe_values (fe_collection[1], quadrature_formula,
+                                     update_values    |  update_gradients |
+                                     update_quadrature_points  |  update_JxW_values);
 
-      if (cell->active_fe_index() == 0)
-       {
-         plain_fe_values.reinit (cell);
+    const unsigned int   n_q_points    = quadrature_formula.size();
 
-         coefficient_values.resize (plain_fe_values.n_quadrature_points);
-         coefficient.value_list (plain_fe_values.get_quadrature_points(),
-                                 coefficient_values);
+    FullMatrix<double>   cell_matrix;
+    Vector<double>       cell_rhs;
 
-         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             {
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
-                 cell_matrix(i,j) += (coefficient_values[q_point] *
-                                      plain_fe_values.shape_grad(i,q_point) *
-                                      plain_fe_values.shape_grad(j,q_point) *
-                                      plain_fe_values.JxW(q_point));
+    std::vector<unsigned int> local_dof_indices;
 
+    const Coefficient<dim> coefficient;
+    std::vector<double>    coefficient_values (n_q_points);
 
-               cell_rhs(i) += (plain_fe_values.shape_value(i,q_point) *
-                               1.0 *
-                               plain_fe_values.JxW(q_point));
-             }
-       }
-      else
-       {
-//TODO: verify that the order of support points equals the order of vertices of the cells, as we use below
-//TODO: remove update_support_points and friends, since they aren't implemented anyway
-         Assert (cell->active_fe_index() == 1, ExcInternalError());
-         Assert (interface_intersects_cell(cell) == true, ExcInternalError());
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
 
-         std::vector<double> level_set_values (GeometryInfo<dim>::vertices_per_cell);
-         for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-           level_set_values[v] = level_set (cell->vertex(v));
+    for (; cell!=endc; ++cell)
+      {
+       const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+       cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+       cell_rhs.reinit (dofs_per_cell);
 
-         FEValues<dim> this_fe_values (fe_collection[1],
-                                       compute_quadrature(quadrature_formula, cell,
-                                                          level_set_values).second,
-                                       update_values    |  update_gradients |
-                                       update_quadrature_points  |  update_JxW_values );
+       cell_matrix = 0;
+       cell_rhs = 0;
 
-         this_fe_values.reinit (cell);
+       if (cell->active_fe_index() == 0)
+         {
+           plain_fe_values.reinit (cell);
 
-         coefficient_values.resize (this_fe_values.n_quadrature_points);
-         coefficient.value_list (this_fe_values.get_quadrature_points(),
-                                 coefficient_values);
+           coefficient_values.resize (plain_fe_values.n_quadrature_points);
+           coefficient.value_list (plain_fe_values.get_quadrature_points(),
+                                   coefficient_values);
 
-         for (unsigned int q_point=0; q_point<this_fe_values.n_quadrature_points; ++q_point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             if (cell->get_fe().system_to_component_index(i).first == 0)
-               {
-                 for (unsigned int j=0; j<dofs_per_cell; ++j)
-                   if (cell->get_fe().system_to_component_index(j).first == 0)
-                     cell_matrix(i,j) += (coefficient_values[q_point] *
-                                          this_fe_values.shape_grad(i,q_point) *
-                                          this_fe_values.shape_grad(j,q_point) *
-                                          this_fe_values.JxW(q_point));
-                   else
-                     cell_matrix(i,j) += (coefficient_values[q_point] *
-                                          this_fe_values.shape_grad(i,q_point)
-                                          *
-                                          ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
-                                            -
-                                            std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))*
-                                           this_fe_values.shape_grad(j,q_point)
-                                           +
-                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
-                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
-                                           this_fe_values.shape_value(j,q_point)) *
-                                          this_fe_values.JxW(q_point));
-
-                 cell_rhs(i) += (this_fe_values.shape_value(i,q_point) *
-                                 1.0 *
-                                 this_fe_values.JxW(q_point));
-               }
-             else
+           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
                {
                  for (unsigned int j=0; j<dofs_per_cell; ++j)
-                   if (cell->get_fe().system_to_component_index(j).first == 0)
-                     cell_matrix(i,j) += (coefficient_values[q_point] *
-                                          ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
-                                            -
-                                            std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
-                                           this_fe_values.shape_grad(i,q_point)
-                                           +
-                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
-                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
-                                           this_fe_values.shape_value(i,q_point)) *
-                                          this_fe_values.shape_grad(j,q_point) *
-                                          this_fe_values.JxW(q_point));
-                   else
-                     cell_matrix(i,j) += (coefficient_values[q_point] *
-                                          ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
-                                            -
-                                            std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
-                                           this_fe_values.shape_grad(i,q_point)
-                                           +
-                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
-                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
-                                           this_fe_values.shape_value(i,q_point)) *
-                                          ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
-                                            -
-                                            std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))*
-                                           this_fe_values.shape_grad(j,q_point)
-                                           +
-                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
-                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
-                                           this_fe_values.shape_value(j,q_point)) *
-                                          this_fe_values.JxW(q_point));
-
-                 cell_rhs(i) += ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
-                                  -
-                                  std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
-                                 this_fe_values.shape_value(i,q_point) *
+                   cell_matrix(i,j) += (coefficient_values[q_point] *
+                                        plain_fe_values.shape_grad(i,q_point) *
+                                        plain_fe_values.shape_grad(j,q_point) *
+                                        plain_fe_values.JxW(q_point));
+
+
+                 cell_rhs(i) += (plain_fe_values.shape_value(i,q_point) *
                                  1.0 *
-                                 this_fe_values.JxW(q_point));
+                                 plain_fe_values.JxW(q_point));
                }
-       }
+         }
+       else
+         {
+//TODO: verify that the order of support points equals the order of vertices of the cells, as we use below
+//TODO: remove update_support_points and friends, since they aren't implemented anyway
+           Assert (cell->active_fe_index() == 1, ExcInternalError());
+           Assert (interface_intersects_cell(cell) == true, ExcInternalError());
+
+           std::vector<double> level_set_values (GeometryInfo<dim>::vertices_per_cell);
+           for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+             level_set_values[v] = level_set (cell->vertex(v));
+
+           FEValues<dim> this_fe_values (fe_collection[1],
+                                         compute_quadrature(quadrature_formula, cell,
+                                                            level_set_values).second,
+                                         update_values    |  update_gradients |
+                                         update_quadrature_points  |  update_JxW_values );
+
+           this_fe_values.reinit (cell);
+
+           coefficient_values.resize (this_fe_values.n_quadrature_points);
+           coefficient.value_list (this_fe_values.get_quadrature_points(),
+                                   coefficient_values);
+
+           for (unsigned int q_point=0; q_point<this_fe_values.n_quadrature_points; ++q_point)
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               if (cell->get_fe().system_to_component_index(i).first == 0)
+                 {
+                   for (unsigned int j=0; j<dofs_per_cell; ++j)
+                     if (cell->get_fe().system_to_component_index(j).first == 0)
+                       cell_matrix(i,j) += (coefficient_values[q_point] *
+                                            this_fe_values.shape_grad(i,q_point) *
+                                            this_fe_values.shape_grad(j,q_point) *
+                                            this_fe_values.JxW(q_point));
+                     else
+                       cell_matrix(i,j) += (coefficient_values[q_point] *
+                                            this_fe_values.shape_grad(i,q_point)
+                                            *
+                                            ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
+                                              -
+                                              std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))*
+                                             this_fe_values.shape_grad(j,q_point)
+                                             +
+                                             grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                             sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                             this_fe_values.shape_value(j,q_point)) *
+                                            this_fe_values.JxW(q_point));
+
+                   cell_rhs(i) += (this_fe_values.shape_value(i,q_point) *
+                                   1.0 *
+                                   this_fe_values.JxW(q_point));
+                 }
+               else
+                 {
+                   for (unsigned int j=0; j<dofs_per_cell; ++j)
+                     if (cell->get_fe().system_to_component_index(j).first == 0)
+                       cell_matrix(i,j) += (coefficient_values[q_point] *
+                                            ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
+                                              -
+                                              std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
+                                             this_fe_values.shape_grad(i,q_point)
+                                             +
+                                             grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                             sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                             this_fe_values.shape_value(i,q_point)) *
+                                            this_fe_values.shape_grad(j,q_point) *
+                                            this_fe_values.JxW(q_point));
+                     else
+                       cell_matrix(i,j) += (coefficient_values[q_point] *
+                                            ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
+                                              -
+                                              std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
+                                             this_fe_values.shape_grad(i,q_point)
+                                             +
+                                             grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                             sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                             this_fe_values.shape_value(i,q_point)) *
+                                            ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
+                                              -
+                                              std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))*
+                                             this_fe_values.shape_grad(j,q_point)
+                                             +
+                                             grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                             sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                             this_fe_values.shape_value(j,q_point)) *
+                                            this_fe_values.JxW(q_point));
+
+                   cell_rhs(i) += ((std::fabs(level_set(this_fe_values.quadrature_point(q_point)))
+                                    -
+                                    std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))*
+                                   this_fe_values.shape_value(i,q_point) *
+                                   1.0 *
+                                   this_fe_values.JxW(q_point));
+                 }
+         }
 
-      local_dof_indices.resize (dofs_per_cell);
-      cell->get_dof_indices (local_dof_indices);
-      constraints.distribute_local_to_global (cell_matrix, cell_rhs,
-                                             local_dof_indices,
-                                             system_matrix, system_rhs);
-    }
+       local_dof_indices.resize (dofs_per_cell);
+       cell->get_dof_indices (local_dof_indices);
+       constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                               local_dof_indices,
+                                               system_matrix, system_rhs);
+      }
 
 
-  std::map<unsigned int,double> boundary_values;
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(2),
-                                           boundary_values);
-  MatrixTools::apply_boundary_values (boundary_values,
-                                     system_matrix,
-                                     solution,
-                                     system_rhs);
+    std::map<unsigned int,double> boundary_values;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             ZeroFunction<dim>(2),
+                                             boundary_values);
+    MatrixTools::apply_boundary_values (boundary_values,
+                                       system_matrix,
+                                       solution,
+                                       system_rhs);
 
-}
+  }
 
 // To integrate the enriched elements we have to find the geometrical decomposition
 // of the original element in subelements. The subelements are used to integrate
@@ -447,654 +449,655 @@ void LaplaceProblem<dim>::assemble_system ()
 // are considered.
 // Type 1: there is not cut. Type 2: a corner of the element is cut. Type 3: two corners are cut.
 
-template <int dim>
-std::pair<unsigned int, Quadrature<dim> >
-LaplaceProblem<dim>::compute_quadrature (const Quadrature<dim> &plain_quadrature,
-                                        const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
-                                        const std::vector<double> &level_set_values                    )
-{
-
-  unsigned int type = 0;
-
-                                  // find the type of cut
-  int sign_ls[GeometryInfo<dim>::vertices_per_cell];
-  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-    {
-      if (level_set_values[v] > 0) sign_ls[v] = 1;
-      else if (level_set_values[v] < 0) sign_ls[v] = -1;
-      else sign_ls[v] = 0;
-    }
+  template <int dim>
+  std::pair<unsigned int, Quadrature<dim> >
+  LaplaceProblem<dim>::compute_quadrature (const Quadrature<dim> &plain_quadrature,
+                                          const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                                          const std::vector<double> &level_set_values                    )
+  {
 
-                                  // the sign of the level set function at the 4 nodes of the elements can be positive + or negative -
-                                  // depending on the sign of the level set function we have the folloing three classes of decomposition
-                                  // type 1: ++++, ----
-                                  // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+
-                                  // type 3: +--+, ++--, +-+-, -++-, --++, -+-+
+    unsigned int type = 0;
 
-  if ( sign_ls[0]==sign_ls[1] & sign_ls[0]==sign_ls[2] & sign_ls[0]==sign_ls[3] ) type =1;
-  else if ( sign_ls[0]*sign_ls[1]*sign_ls[2]*sign_ls[3] < 0 ) type = 2;
-  else type = 3;
+                                    // find the type of cut
+    int sign_ls[GeometryInfo<dim>::vertices_per_cell];
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+      {
+       if (level_set_values[v] > 0) sign_ls[v] = 1;
+       else if (level_set_values[v] < 0) sign_ls[v] = -1;
+       else sign_ls[v] = 0;
+      }
 
-  unsigned int Pos = 100;
+                                    // the sign of the level set function at the 4 nodes of the elements can be positive + or negative -
+                                    // depending on the sign of the level set function we have the folloing three classes of decomposition
+                                    // type 1: ++++, ----
+                                    // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+
+                                    // type 3: +--+, ++--, +-+-, -++-, --++, -+-+
 
-  Point<dim> v0(0,0);
-  Point<dim> v1(1,0);
-  Point<dim> v2(0,1);
-  Point<dim> v3(1,1);
+    if ( sign_ls[0]==sign_ls[1] & sign_ls[0]==sign_ls[2] & sign_ls[0]==sign_ls[3] ) type =1;
+    else if ( sign_ls[0]*sign_ls[1]*sign_ls[2]*sign_ls[3] < 0 ) type = 2;
+    else type = 3;
 
-  Point<dim> A(0,0);
-  Point<dim> B(0,0);
-  Point<dim> C(0,0);
-  Point<dim> D(0,0);
-  Point<dim> E(0,0);
-  Point<dim> F(0,0);
-
-  if (type == 1)
-    return std::pair<unsigned int, Quadrature<dim> >(1, plain_quadrature);
-
-  if (type==2)
-    {
-      const unsigned int   n_q_points    = plain_quadrature.size();
-
-                                      // loop over all subelements for integration
-                                      // in type 2 there are 5 subelements
-
-      Quadrature<dim> xfem_quadrature(5*n_q_points);
-
-      std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
-
-      if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0;
-      else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1;
-      else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2;
-      else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3;
-      else assert(0); // error message
-
-                                      // Find cut coordinates
-
-                                      // deal.ii local coordinates
-
-                                      //    2-------3
-                                      //    |       |
-                                      //               |       |
-                                      //               |       |
-                                      //               0-------1
-
-      if (Pos == 0)
-       {
-         A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]);
-         B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]);
-         A(1) = 0.;
-         B(0) = 0.;
-         C(0) = 0.5*( A(0) + B(0) );
-         C(1) = 0.5*( A(1) + B(1) );
-         D(0) = 2./3. * C(0);
-         D(1) = 2./3. * C(1);
-         E(0) = 0.5*A(0);
-         E(1) = 0.;
-         F(0) = 0.;
-         F(1) = 0.5*B(1);
-       }
-      else if (Pos == 1)
-       {
-         A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]);
-         B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]);
-         A(1) = 0.;
-         B(0) = 1.;
-         C(0) = 0.5*( A(0) + B(0) );
-         C(1) = 0.5*( A(1) + B(1) );
-         D(0) = 1./3. + 2./3. * C(0);
-         D(1) = 2./3. * C(1);
-         E(0) = 0.5*(1 + A(0));
-         E(1) = 0.;
-         F(0) = 1.;
-         F(1) = 0.5*B(1);
-       }
-      else if (Pos == 2)
-       {
-         A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]);
-         B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]);
-         A(1) = 1.;
-         B(0) = 0.;
-         C(0) = 0.5*( A(0) + B(0) );
-         C(1) = 0.5*( A(1) + B(1) );
-         D(0) = 2./3. * C(0);
-         D(1) = 1./3. + 2./3. * C(1);
-         E(0) = 0.5* A(0);
-         E(1) = 1.;
-         F(0) = 0.;
-         F(1) = 0.5*( 1. + B(1) );
-       }
-      else if (Pos == 3)
-       {
-         A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]);
-         B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]);
-         A(1) = 1.;
-         B(0) = 1.;
-         C(0) = 0.5*( A(0) + B(0) );
-         C(1) = 0.5*( A(1) + B(1) );
-         D(0) = 1./3. + 2./3. * C(0);
-         D(1) = 1./3. + 2./3. * C(1);
-         E(0) = 0.5*( 1. + A(0) );
-         E(1) = 1.;
-         F(0) = 1.;
-         F(1) = 0.5*( 1. + B(1) );
-       }
-
-                                      //std::cout << A << std::endl;
-                                      //std::cout << B << std::endl;
-                                      //std::cout << C << std::endl;
-                                      //std::cout << D << std::endl;
-                                      //std::cout << E << std::endl;
-                                      //std::cout << F << std::endl;
-
-      std::string filename = "vertices.dat";
-      std::ofstream output (filename.c_str());
-      output << "#vertices of xfem subcells" << std::endl;
-      output << v0(0) << "   " << v0(1) << std::endl;
-      output << v1(0) << "   " << v1(1) << std::endl;
-      output << v3(0) << "   " << v3(1) << std::endl;
-      output << v2(0) << "   " << v2(1) << std::endl;
-      output << std::endl;
-      output << A(0) << "   " << A(1) << std::endl;
-      output << B(0) << "   " << B(1) << std::endl;
-      output << std::endl;
-      output << C(0) << "   " << C(1) << std::endl;
-      output << D(0) << "   " << D(1) << std::endl;
-      output << std::endl;
-      output << D(0) << "   " << D(1) << std::endl;
-      output << E(0) << "   " << E(1) << std::endl;
-      output << std::endl;
-      output << D(0) << "   " << D(1) << std::endl;
-      output << F(0) << "   " << F(1) << std::endl;
-      output << std::endl;
-
-      if (Pos==0)
-       output << v3(0) << "   " << v3(1) << std::endl;
-      else if (Pos==1)
-       output << v2(0) << "   " << v2(1) << std::endl;
-      else if (Pos==2)
-       output << v1(0) << "   " << v1(1) << std::endl;
-      else if (Pos==3)
-       output << v0(0) << "   " << v0(1) << std::endl;
-      output << C(0) << "   " << C(1) << std::endl;
-
-      Point<dim> subcell_vertices[10];
-      subcell_vertices[0] = v0;
-      subcell_vertices[1] = v1;
-      subcell_vertices[2] = v2;
-      subcell_vertices[3] = v3;
-      subcell_vertices[4] = A;
-      subcell_vertices[5] = B;
-      subcell_vertices[6] = C;
-      subcell_vertices[7] = D;
-      subcell_vertices[8] = E;
-      subcell_vertices[9] = F;
-
-      std::vector<Point<dim> > xfem_points;
-      std::vector<double>      xfem_weights;
-
-                                      // lookup table for the decomposition
-
-      if (dim==2)
-       {
-         unsigned int subcell_v_indices[4][5][4] = {
-               {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}},
-               {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}},
-               {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}},
-               {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}}
-         };
-
-         for (unsigned int subcell = 0; subcell<5; subcell++)
-           {
-                                              //std::cout << "subcell   : " << subcell << std::endl;
-             std::vector<Point<dim> > vertices;
-             for (unsigned int i=0; i<4; i++)
-               {
-                 vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
-                                                  //std::cout << "i         : " << i << std::endl;
-                                                  //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
-                                                  //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
-               }
-                                              //std::cout << std::endl;
-                                              // create quadrature rule
-             append_quadrature( plain_quadrature,
-                                vertices,
-                                xfem_points,
-                                xfem_weights);
-                                              //initialize xfem_quadrature with quadrature points of all subelements
-             xfem_quadrature.initialize(xfem_points, xfem_weights);
-           }
-       }
+    unsigned int Pos = 100;
 
-      Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError());
-      return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
-    }
+    Point<dim> v0(0,0);
+    Point<dim> v1(1,0);
+    Point<dim> v2(0,1);
+    Point<dim> v3(1,1);
 
-                                  // Type three decomposition
-                                  // (+--+, ++--, +-+-, -++-, --++, -+-+)
+    Point<dim> A(0,0);
+    Point<dim> B(0,0);
+    Point<dim> C(0,0);
+    Point<dim> D(0,0);
+    Point<dim> E(0,0);
+    Point<dim> F(0,0);
 
-  if (type==3)
-    {
-      const unsigned int   n_q_points    = plain_quadrature.size();
-
-                                      // loop over all subelements for integration
-                                      // in type 2 there are 5 subelements
-
-      Quadrature<dim> xfem_quadrature(5*n_q_points);
-
-      std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
-
-      if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] )
-       {
-         Pos = 0;
-         A(0) = 0.;
-         A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2]));
-         B(0) = 1.;
-         B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3]));
-       }
-      else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] )
-       {
-         Pos = 1;
-         A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1]));
-         A(1) = 0.;
-         B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3]));
-         B(1) = 1.;
-       }
-      else if ( sign_ls[0]==sign_ls[3] && sign_ls[1]==sign_ls[2] )
-       {
-         std::cout << "Error: the element has two cut lines and this is not allowed" << std::endl;
-         assert(0);
-       }
-      else
-       {
-         std::cout << "Error: the level set function has not the right values" << std::endl;
-         assert(0);
-       }
-
-                                      //std::cout << "Pos " << Pos << std::endl;
-                                      //std::cout << A << std::endl;
-                                      //std::cout << B << std::endl;
-      std::string filename = "vertices.dat";
-      std::ofstream output (filename.c_str());
-      output << "#vertices of xfem subcells" << std::endl;
-      output << A(0) << "   " << A(1) << std::endl;
-      output << B(0) << "   " << B(1) << std::endl;
-
-                                      //fill xfem_quadrature
-      Point<dim> subcell_vertices[6];
-      subcell_vertices[0] = v0;
-      subcell_vertices[1] = v1;
-      subcell_vertices[2] = v2;
-      subcell_vertices[3] = v3;
-      subcell_vertices[4] = A;
-      subcell_vertices[5] = B;
-
-      std::vector<Point<dim> > xfem_points;
-      std::vector<double>      xfem_weights;
-
-      if (dim==2)
-       {
-         unsigned int subcell_v_indices[2][2][4] = {
-               {{0,1,4,5}, {4,5,2,3}},
-               {{0,4,2,5}, {4,1,5,3}}
-         };
-
-                                          //std::cout << "Pos       : " << Pos << std::endl;
-         for (unsigned int subcell = 0; subcell<2; subcell++)
-           {
-                                              //std::cout << "subcell   : " << subcell << std::endl;
-             std::vector<Point<dim> > vertices;
-             for (unsigned int i=0; i<4; i++)
-               {
-                 vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
-                                                  //std::cout << "i         : " << i << std::endl;
-                                                  //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
-                                                  //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
-               }
-                                              //std::cout << std::endl;
-                                              // create quadrature rule
-             append_quadrature( plain_quadrature,
-                                vertices,
-                                xfem_points,
-                                xfem_weights);
-                                              //initialize xfem_quadrature with quadrature points of all subelements
-             xfem_quadrature.initialize(xfem_points, xfem_weights);
-           }
-       }
-      Assert (xfem_quadrature.size() == plain_quadrature.size() * 2, ExcInternalError());
-      return std::pair<unsigned int, Quadrature<dim> >(3, xfem_quadrature);
-    }
+    if (type == 1)
+      return std::pair<unsigned int, Quadrature<dim> >(1, plain_quadrature);
 
-  return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
+    if (type==2)
+      {
+       const unsigned int   n_q_points    = plain_quadrature.size();
 
-}
+                                        // loop over all subelements for integration
+                                        // in type 2 there are 5 subelements
 
-template <int dim>
-void LaplaceProblem<dim>::append_quadrature ( const Quadrature<dim> &plain_quadrature,
-                                             const std::vector<Point<dim> > &v,
-                                             std::vector<Point<dim> > &xfem_points,
-                                             std::vector<double>      &xfem_weights)
+       Quadrature<dim> xfem_quadrature(5*n_q_points);
 
-{
-                                  // Project integration points into sub-elements.
-                                  // This maps quadrature points from a reference element to a subelement of a reference element.
-                                  // To implement the action of this map the coordinates of the subelements have been calculated (A(0)...F(0),A(1)...F(1))
-                                  // the coordinates of the quadrature points are given by the bi-linear map defined by the form functions
-                                  // $x^\prime_i = \sum_j v^\prime \phi_j(x^hat_i)$, where the $\phi_j$ are the shape functions of the FEQ.
+       std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
 
-  unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
+       if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0;
+       else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1;
+       else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2;
+       else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3;
+       else assert(0); // error message
 
-  std::vector<Point<dim> > q_points = plain_quadrature.get_points();
-  std::vector<Point<dim> > q_transf(q_points.size());
-  std::vector<double> W = plain_quadrature.get_weights();
-  std::vector<double> phi(n_v);
-  std::vector<Tensor<1,dim> > grad_phi(n_v);
+                                        // Find cut coordinates
 
-  const unsigned int   n_q_points    = plain_quadrature.size();
+                                        // deal.ii local coordinates
 
-  std::vector<double> JxW(n_q_points);
+                                        //    2-------3
+                                        //    |       |
+                                        //             |       |
+                                        //             |       |
+                                        //             0-------1
 
-  for ( unsigned int i = 0; i < n_q_points; i++)
-    {
-      switch (dim)
-       {
-         case 2:
+       if (Pos == 0)
          {
-           double xi  = q_points[i](0);
-           double eta = q_points[i](1);
-
-                                            // Define shape functions on reference element
-                                            // we consider a bi-linear mapping
-           phi[0] = (1. - xi) * (1. - eta);
-           phi[1] = xi * (1. - eta);
-           phi[2] = (1. - xi) * eta;
-           phi[3] = xi * eta;
-
-           grad_phi[0][0] = (-1. + eta);
-           grad_phi[1][0] = (1. - eta);
-           grad_phi[2][0] = -eta;
-           grad_phi[3][0] = eta;
-
-           grad_phi[0][1] = (-1. + xi);
-           grad_phi[1][1] = -xi;
-           grad_phi[2][1] = 1-xi;
-           grad_phi[3][1] = xi;
-
-           break;
+           A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]);
+           B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]);
+           A(1) = 0.;
+           B(0) = 0.;
+           C(0) = 0.5*( A(0) + B(0) );
+           C(1) = 0.5*( A(1) + B(1) );
+           D(0) = 2./3. * C(0);
+           D(1) = 2./3. * C(1);
+           E(0) = 0.5*A(0);
+           E(1) = 0.;
+           F(0) = 0.;
+           F(1) = 0.5*B(1);
+         }
+       else if (Pos == 1)
+         {
+           A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]);
+           B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]);
+           A(1) = 0.;
+           B(0) = 1.;
+           C(0) = 0.5*( A(0) + B(0) );
+           C(1) = 0.5*( A(1) + B(1) );
+           D(0) = 1./3. + 2./3. * C(0);
+           D(1) = 2./3. * C(1);
+           E(0) = 0.5*(1 + A(0));
+           E(1) = 0.;
+           F(0) = 1.;
+           F(1) = 0.5*B(1);
+         }
+       else if (Pos == 2)
+         {
+           A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]);
+           B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]);
+           A(1) = 1.;
+           B(0) = 0.;
+           C(0) = 0.5*( A(0) + B(0) );
+           C(1) = 0.5*( A(1) + B(1) );
+           D(0) = 2./3. * C(0);
+           D(1) = 1./3. + 2./3. * C(1);
+           E(0) = 0.5* A(0);
+           E(1) = 1.;
+           F(0) = 0.;
+           F(1) = 0.5*( 1. + B(1) );
+         }
+       else if (Pos == 3)
+         {
+           A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]);
+           B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]);
+           A(1) = 1.;
+           B(0) = 1.;
+           C(0) = 0.5*( A(0) + B(0) );
+           C(1) = 0.5*( A(1) + B(1) );
+           D(0) = 1./3. + 2./3. * C(0);
+           D(1) = 1./3. + 2./3. * C(1);
+           E(0) = 0.5*( 1. + A(0) );
+           E(1) = 1.;
+           F(0) = 1.;
+           F(1) = 0.5*( 1. + B(1) );
          }
 
-         default:
-               Assert (false, ExcNotImplemented());
-       }
-
-
-      Tensor<2,dim> jacobian;
+                                        //std::cout << A << std::endl;
+                                        //std::cout << B << std::endl;
+                                        //std::cout << C << std::endl;
+                                        //std::cout << D << std::endl;
+                                        //std::cout << E << std::endl;
+                                        //std::cout << F << std::endl;
 
-                                      // Calculate Jacobian of transformation
-      for (unsigned int d=0; d<dim; ++d)
-       for (unsigned int e=0; e<dim; ++e)
+       std::string filename = "vertices.dat";
+       std::ofstream output (filename.c_str());
+       output << "#vertices of xfem subcells" << std::endl;
+       output << v0(0) << "   " << v0(1) << std::endl;
+       output << v1(0) << "   " << v1(1) << std::endl;
+       output << v3(0) << "   " << v3(1) << std::endl;
+       output << v2(0) << "   " << v2(1) << std::endl;
+       output << std::endl;
+       output << A(0) << "   " << A(1) << std::endl;
+       output << B(0) << "   " << B(1) << std::endl;
+       output << std::endl;
+       output << C(0) << "   " << C(1) << std::endl;
+       output << D(0) << "   " << D(1) << std::endl;
+       output << std::endl;
+       output << D(0) << "   " << D(1) << std::endl;
+       output << E(0) << "   " << E(1) << std::endl;
+       output << std::endl;
+       output << D(0) << "   " << D(1) << std::endl;
+       output << F(0) << "   " << F(1) << std::endl;
+       output << std::endl;
+
+       if (Pos==0)
+         output << v3(0) << "   " << v3(1) << std::endl;
+       else if (Pos==1)
+         output << v2(0) << "   " << v2(1) << std::endl;
+       else if (Pos==2)
+         output << v1(0) << "   " << v1(1) << std::endl;
+       else if (Pos==3)
+         output << v0(0) << "   " << v0(1) << std::endl;
+       output << C(0) << "   " << C(1) << std::endl;
+
+       Point<dim> subcell_vertices[10];
+       subcell_vertices[0] = v0;
+       subcell_vertices[1] = v1;
+       subcell_vertices[2] = v2;
+       subcell_vertices[3] = v3;
+       subcell_vertices[4] = A;
+       subcell_vertices[5] = B;
+       subcell_vertices[6] = C;
+       subcell_vertices[7] = D;
+       subcell_vertices[8] = E;
+       subcell_vertices[9] = F;
+
+       std::vector<Point<dim> > xfem_points;
+       std::vector<double>      xfem_weights;
+
+                                        // lookup table for the decomposition
+
+       if (dim==2)
          {
-           for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
+           unsigned int subcell_v_indices[4][5][4] = {
+                 {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}},
+                 {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}},
+                 {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}},
+                 {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}}
+           };
+
+           for (unsigned int subcell = 0; subcell<5; subcell++)
              {
-               jacobian[d][e] += grad_phi[j][e] * v[j](d);
+                                                //std::cout << "subcell   : " << subcell << std::endl;
+               std::vector<Point<dim> > vertices;
+               for (unsigned int i=0; i<4; i++)
+                 {
+                   vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
+                                                    //std::cout << "i         : " << i << std::endl;
+                                                    //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
+                                                    //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
+                 }
+                                                //std::cout << std::endl;
+                                                // create quadrature rule
+               append_quadrature( plain_quadrature,
+                                  vertices,
+                                  xfem_points,
+                                  xfem_weights);
+                                                //initialize xfem_quadrature with quadrature points of all subelements
+               xfem_quadrature.initialize(xfem_points, xfem_weights);
              }
          }
 
-      double detJ = determinant(jacobian);
-      xfem_weights.push_back (W[i] * detJ);
+       Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError());
+       return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
+      }
 
-                                      // Map integration points from reference element to subcell of reference element
-      Point<dim> q_prime;
-      for (unsigned int d=0; d<dim; ++d)
-       for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
-         q_prime[d] += v[j](d) * phi[j];
-      xfem_points.push_back(q_prime);
-    }
+                                    // Type three decomposition
+                                    // (+--+, ++--, +-+-, -++-, --++, -+-+)
 
-}
+    if (type==3)
+      {
+       const unsigned int   n_q_points    = plain_quadrature.size();
 
+                                        // loop over all subelements for integration
+                                        // in type 2 there are 5 subelements
 
-template <int dim>
-void LaplaceProblem<dim>::solve ()
-{
-  SolverControl           solver_control (1000, 1e-12);
-  SolverCG<>              solver (solver_control);
+       Quadrature<dim> xfem_quadrature(5*n_q_points);
+
+       std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
 
-  PreconditionSSOR<> preconditioner;
-  preconditioner.initialize(system_matrix, 1.2);
+       if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] )
+         {
+           Pos = 0;
+           A(0) = 0.;
+           A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2]));
+           B(0) = 1.;
+           B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3]));
+         }
+       else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] )
+         {
+           Pos = 1;
+           A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1]));
+           A(1) = 0.;
+           B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3]));
+           B(1) = 1.;
+         }
+       else if ( sign_ls[0]==sign_ls[3] && sign_ls[1]==sign_ls[2] )
+         {
+           std::cout << "Error: the element has two cut lines and this is not allowed" << std::endl;
+           assert(0);
+         }
+       else
+         {
+           std::cout << "Error: the level set function has not the right values" << std::endl;
+           assert(0);
+         }
 
-  solver.solve (system_matrix, solution, system_rhs,
-               preconditioner);
+                                        //std::cout << "Pos " << Pos << std::endl;
+                                        //std::cout << A << std::endl;
+                                        //std::cout << B << std::endl;
+       std::string filename = "vertices.dat";
+       std::ofstream output (filename.c_str());
+       output << "#vertices of xfem subcells" << std::endl;
+       output << A(0) << "   " << A(1) << std::endl;
+       output << B(0) << "   " << B(1) << std::endl;
+
+                                        //fill xfem_quadrature
+       Point<dim> subcell_vertices[6];
+       subcell_vertices[0] = v0;
+       subcell_vertices[1] = v1;
+       subcell_vertices[2] = v2;
+       subcell_vertices[3] = v3;
+       subcell_vertices[4] = A;
+       subcell_vertices[5] = B;
+
+       std::vector<Point<dim> > xfem_points;
+       std::vector<double>      xfem_weights;
+
+       if (dim==2)
+         {
+           unsigned int subcell_v_indices[2][2][4] = {
+                 {{0,1,4,5}, {4,5,2,3}},
+                 {{0,4,2,5}, {4,1,5,3}}
+           };
 
-  constraints.distribute (solution);
-}
+                                            //std::cout << "Pos       : " << Pos << std::endl;
+           for (unsigned int subcell = 0; subcell<2; subcell++)
+             {
+                                                //std::cout << "subcell   : " << subcell << std::endl;
+               std::vector<Point<dim> > vertices;
+               for (unsigned int i=0; i<4; i++)
+                 {
+                   vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
+                                                    //std::cout << "i         : " << i << std::endl;
+                                                    //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
+                                                    //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
+                 }
+                                                //std::cout << std::endl;
+                                                // create quadrature rule
+               append_quadrature( plain_quadrature,
+                                  vertices,
+                                  xfem_points,
+                                  xfem_weights);
+                                                //initialize xfem_quadrature with quadrature points of all subelements
+               xfem_quadrature.initialize(xfem_points, xfem_weights);
+             }
+         }
+       Assert (xfem_quadrature.size() == plain_quadrature.size() * 2, ExcInternalError());
+       return std::pair<unsigned int, Quadrature<dim> >(3, xfem_quadrature);
+      }
 
+    return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
 
+  }
 
-template <int dim>
-void LaplaceProblem<dim>::refine_grid ()
-{
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  template <int dim>
+  void LaplaceProblem<dim>::append_quadrature ( const Quadrature<dim> &plain_quadrature,
+                                               const std::vector<Point<dim> > &v,
+                                               std::vector<Point<dim> > &xfem_points,
+                                               std::vector<double>      &xfem_weights)
 
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                     QGauss<dim-1>(3),
-                                     typename FunctionMap<dim>::type(),
-                                     solution,
-                                     estimated_error_per_cell);
+  {
+                                    // Project integration points into sub-elements.
+                                    // This maps quadrature points from a reference element to a subelement of a reference element.
+                                    // To implement the action of this map the coordinates of the subelements have been calculated (A(0)...F(0),A(1)...F(1))
+                                    // the coordinates of the quadrature points are given by the bi-linear map defined by the form functions
+                                    // $x^\prime_i = \sum_j v^\prime \phi_j(x^hat_i)$, where the $\phi_j$ are the shape functions of the FEQ.
 
-  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                  estimated_error_per_cell,
-                                                  0.3, 0.03);
+    unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
 
-  triangulation.execute_coarsening_and_refinement ();
-}
+    std::vector<Point<dim> > q_points = plain_quadrature.get_points();
+    std::vector<Point<dim> > q_transf(q_points.size());
+    std::vector<double> W = plain_quadrature.get_weights();
+    std::vector<double> phi(n_v);
+    std::vector<Tensor<1,dim> > grad_phi(n_v);
 
+    const unsigned int   n_q_points    = plain_quadrature.size();
 
+    std::vector<double> JxW(n_q_points);
 
-template <int dim>
-class Postprocessor : public DataPostprocessor<dim>
-{
-  public:
-    virtual
-    void
-    compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
-                                      const std::vector<std::vector<Tensor<1,dim> > > &duh,
-                                      const std::vector<std::vector<Tensor<2,dim> > > &dduh,
-                                      const std::vector<Point<dim> >                  &normals,
-                                      const std::vector<Point<dim> >                  &evaluation_points,
-                                      std::vector<Vector<double> >                    &computed_quantities) const;
+    for ( unsigned int i = 0; i < n_q_points; i++)
+      {
+       switch (dim)
+         {
+           case 2:
+           {
+             double xi  = q_points[i](0);
+             double eta = q_points[i](1);
+
+                                              // Define shape functions on reference element
+                                              // we consider a bi-linear mapping
+             phi[0] = (1. - xi) * (1. - eta);
+             phi[1] = xi * (1. - eta);
+             phi[2] = (1. - xi) * eta;
+             phi[3] = xi * eta;
+
+             grad_phi[0][0] = (-1. + eta);
+             grad_phi[1][0] = (1. - eta);
+             grad_phi[2][0] = -eta;
+             grad_phi[3][0] = eta;
+
+             grad_phi[0][1] = (-1. + xi);
+             grad_phi[1][1] = -xi;
+             grad_phi[2][1] = 1-xi;
+             grad_phi[3][1] = xi;
+
+             break;
+           }
 
-    virtual std::vector<std::string> get_names () const;
+           default:
+                 Assert (false, ExcNotImplemented());
+         }
 
-    virtual unsigned int n_output_variables() const;
 
-    virtual
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    get_data_component_interpretation () const;
+       Tensor<2,dim> jacobian;
 
-    virtual UpdateFlags get_needed_update_flags () const;
-};
+                                        // Calculate Jacobian of transformation
+       for (unsigned int d=0; d<dim; ++d)
+         for (unsigned int e=0; e<dim; ++e)
+           {
+             for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
+               {
+                 jacobian[d][e] += grad_phi[j][e] * v[j](d);
+               }
+           }
 
+       double detJ = determinant(jacobian);
+       xfem_weights.push_back (W[i] * detJ);
 
-template <int dim>
-std::vector<std::string>
-Postprocessor<dim>::get_names() const
-{
-  std::vector<std::string> solution_names (1, "total_solution");
-  solution_names.push_back ("error");
-  return solution_names;
-}
+                                        // Map integration points from reference element to subcell of reference element
+       Point<dim> q_prime;
+       for (unsigned int d=0; d<dim; ++d)
+         for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
+           q_prime[d] += v[j](d) * phi[j];
+       xfem_points.push_back(q_prime);
+      }
 
+  }
 
-template <int dim>
-unsigned int
-Postprocessor<dim>::n_output_variables() const
-{
-  return get_names().size();
-}
 
+  template <int dim>
+  void LaplaceProblem<dim>::solve ()
+  {
+    SolverControl           solver_control (1000, 1e-12);
+    SolverCG<>              solver (solver_control);
 
-template <int dim>
-std::vector<DataComponentInterpretation::DataComponentInterpretation>
-Postprocessor<dim>::
-get_data_component_interpretation () const
-{
-  std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    interpretation (2,
-                   DataComponentInterpretation::component_is_scalar);
-  return interpretation;
-}
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
 
+    solver.solve (system_matrix, solution, system_rhs,
+                 preconditioner);
 
-template <int dim>
-UpdateFlags
-Postprocessor<dim>::get_needed_update_flags() const
-{
-  return update_values | update_q_points;
-}
+    constraints.distribute (solution);
+  }
 
 
-template <int dim>
-void
-Postprocessor<dim>::
-compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
-                                  const std::vector<std::vector<Tensor<1,dim> > > &/*duh*/,
-                                  const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
-                                  const std::vector<Point<dim> >                  &/*normals*/,
-                                  const std::vector<Point<dim> >                  &evaluation_points,
-                                  std::vector<Vector<double> >                    &computed_quantities) const
-{
-  const unsigned int n_quadrature_points = uh.size();
-  Assert (computed_quantities.size() == n_quadrature_points,  ExcInternalError());
-  Assert (uh[0].size() == 2,                                  ExcInternalError());
-  Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError());
 
-  for (unsigned int q=0; q<n_quadrature_points; ++q)
-    {
-      computed_quantities[q](0)
-       = (uh[q](0)
-          +
-          uh[q](1) * std::fabs(level_set(evaluation_points[q])));
-      computed_quantities[q](1)
-       = (computed_quantities[q](0)
-          -
-          exact_solution (evaluation_points[q]));
-    }
-}
+  template <int dim>
+  void LaplaceProblem<dim>::refine_grid ()
+  {
+    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                       QGauss<dim-1>(3),
+                                       typename FunctionMap<dim>::type(),
+                                       solution,
+                                       estimated_error_per_cell);
 
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.3, 0.03);
 
-template <int dim>
-void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
-{
-  Assert (cycle < 10, ExcNotImplemented());
+    triangulation.execute_coarsening_and_refinement ();
+  }
 
-  std::string filename = "solution-";
-  filename += ('0' + cycle);
-  //filename += ".vtk";
-  filename += ".gmv";
 
-  std::ofstream output (filename.c_str());
 
-  Postprocessor<dim> postprocessor;
-  DataOut<dim,hp::DoFHandler<dim> > data_out;
+  template <int dim>
+  class Postprocessor : public DataPostprocessor<dim>
+  {
+    public:
+      virtual
+      void
+      compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                        const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                        const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+                                        const std::vector<Point<dim> >                  &normals,
+                                        const std::vector<Point<dim> >                  &evaluation_points,
+                                        std::vector<Vector<double> >                    &computed_quantities) const;
 
-  data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (solution, "solution");
-  data_out.add_data_vector (solution, postprocessor);
-  data_out.build_patches (5);
+      virtual std::vector<std::string> get_names () const;
 
-  //data_out.write_vtk (output);
-  data_out.write_gmv (output);
-}
+      virtual unsigned int n_output_variables() const;
 
+      virtual
+      std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      get_data_component_interpretation () const;
 
+      virtual UpdateFlags get_needed_update_flags () const;
+  };
 
-template <int dim>
-void LaplaceProblem<dim>::compute_error () const
-{
-  hp::QCollection<dim> q_collection;
-  q_collection.push_back (QGauss<dim>(2));
-  q_collection.push_back (QIterated<dim>(QGauss<1>(2), 4));
 
-  hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
-                                 update_values | update_q_points | update_JxW_values);
+  template <int dim>
+  std::vector<std::string>
+  Postprocessor<dim>::get_names() const
+  {
+    std::vector<std::string> solution_names (1, "total_solution");
+    solution_names.push_back ("error");
+    return solution_names;
+  }
 
-  double l2_error_square = 0;
 
-  std::vector<Vector<double> > solution_values;
+  template <int dim>
+  unsigned int
+  Postprocessor<dim>::n_output_variables() const
+  {
+    return get_names().size();
+  }
 
-  typename hp::DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
 
-  for (; cell!=endc; ++cell)
-    {
-      hp_fe_values.reinit (cell);
-
-      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
-
-      solution_values.resize (fe_values.n_quadrature_points,
-                             Vector<double>(2));
-      fe_values.get_function_values (solution,
-                                    solution_values);
-
-      for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
-       {
-         const double local_error = (solution_values[q](0)
-                                     +
-                                     std::fabs(level_set(fe_values.quadrature_point(q))) *
-                                     solution_values[q](1)
-                                     -
-                                     exact_solution (fe_values.quadrature_point(q)));
-         l2_error_square += local_error * local_error * fe_values.JxW(q);
-       }
-    }
+  template <int dim>
+  std::vector<DataComponentInterpretation::DataComponentInterpretation>
+  Postprocessor<dim>::
+  get_data_component_interpretation () const
+  {
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      interpretation (2,
+                     DataComponentInterpretation::component_is_scalar);
+    return interpretation;
+  }
+
+
+  template <int dim>
+  UpdateFlags
+  Postprocessor<dim>::get_needed_update_flags() const
+  {
+    return update_values | update_q_points;
+  }
+
+
+  template <int dim>
+  void
+  Postprocessor<dim>::
+  compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                    const std::vector<std::vector<Tensor<1,dim> > > &/*duh*/,
+                                    const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+                                    const std::vector<Point<dim> >                  &/*normals*/,
+                                    const std::vector<Point<dim> >                  &evaluation_points,
+                                    std::vector<Vector<double> >                    &computed_quantities) const
+  {
+    const unsigned int n_quadrature_points = uh.size();
+    Assert (computed_quantities.size() == n_quadrature_points,  ExcInternalError());
+    Assert (uh[0].size() == 2,                                  ExcInternalError());
+    Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError());
+
+    for (unsigned int q=0; q<n_quadrature_points; ++q)
+      {
+       computed_quantities[q](0)
+         = (uh[q](0)
+            +
+            uh[q](1) * std::fabs(level_set(evaluation_points[q])));
+       computed_quantities[q](1)
+         = (computed_quantities[q](0)
+            -
+            exact_solution (evaluation_points[q]));
+      }
+  }
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+  {
+    Assert (cycle < 10, ExcNotImplemented());
+
+    std::string filename = "solution-";
+    filename += ('0' + cycle);
+                                    //filename += ".vtk";
+    filename += ".gmv";
+
+    std::ofstream output (filename.c_str());
+
+    Postprocessor<dim> postprocessor;
+    DataOut<dim,hp::DoFHandler<dim> > data_out;
+
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, "solution");
+    data_out.add_data_vector (solution, postprocessor);
+    data_out.build_patches (5);
+
+                                    //data_out.write_vtk (output);
+    data_out.write_gmv (output);
+  }
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::compute_error () const
+  {
+    hp::QCollection<dim> q_collection;
+    q_collection.push_back (QGauss<dim>(2));
+    q_collection.push_back (QIterated<dim>(QGauss<1>(2), 4));
+
+    hp::FEValues<dim> hp_fe_values (fe_collection, q_collection,
+                                   update_values | update_q_points | update_JxW_values);
+
+    double l2_error_square = 0;
+
+    std::vector<Vector<double> > solution_values;
+
+    typename hp::DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+
+    for (; cell!=endc; ++cell)
+      {
+       hp_fe_values.reinit (cell);
+
+       const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
+
+       solution_values.resize (fe_values.n_quadrature_points,
+                               Vector<double>(2));
+       fe_values.get_function_values (solution,
+                                      solution_values);
+
+       for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+         {
+           const double local_error = (solution_values[q](0)
+                                       +
+                                       std::fabs(level_set(fe_values.quadrature_point(q))) *
+                                       solution_values[q](1)
+                                       -
+                                       exact_solution (fe_values.quadrature_point(q)));
+           l2_error_square += local_error * local_error * fe_values.JxW(q);
+         }
+      }
 
-  std::cout << "   L2 error = " << std::sqrt (l2_error_square)
-           << std::endl;
-}
+    std::cout << "   L2 error = " << std::sqrt (l2_error_square)
+             << std::endl;
+  }
 
 
 
 
-template <int dim>
-void LaplaceProblem<dim>::run ()
-{
-  for (unsigned int cycle=0; cycle<6; ++cycle)
-    {
-      std::cout << "Cycle " << cycle << ':' << std::endl;
+  template <int dim>
+  void LaplaceProblem<dim>::run ()
+  {
+    for (unsigned int cycle=0; cycle<6; ++cycle)
+      {
+       std::cout << "Cycle " << cycle << ':' << std::endl;
 
-      if (cycle == 0)
-       {
-         GridGenerator::hyper_ball (triangulation);
-               //GridGenerator::hyper_cube (triangulation, -1, 1);
+       if (cycle == 0)
+         {
+           GridGenerator::hyper_ball (triangulation);
+                                            //GridGenerator::hyper_cube (triangulation, -1, 1);
 
-         static const HyperBallBoundary<dim> boundary;
-         triangulation.set_boundary (0, boundary);
+           static const HyperBallBoundary<dim> boundary;
+           triangulation.set_boundary (0, boundary);
 
-         triangulation.refine_global (2);
-       }
-      else
-       triangulation.refine_global (1);
+           triangulation.refine_global (2);
+         }
+       else
+         triangulation.refine_global (1);
 //     refine_grid ();
 
 
-      std::cout << "   Number of active cells:       "
-               << triangulation.n_active_cells()
-               << std::endl;
+       std::cout << "   Number of active cells:       "
+                 << triangulation.n_active_cells()
+                 << std::endl;
 
-      setup_system ();
+       setup_system ();
 
-      std::cout << "   Number of degrees of freedom: "
-               << dof_handler.n_dofs()
-               << std::endl;
+       std::cout << "   Number of degrees of freedom: "
+                 << dof_handler.n_dofs()
+                 << std::endl;
 
-      assemble_system ();
-      solve ();
-      compute_error ();
-      output_results (cycle);
-    }
+       assemble_system ();
+       solve ();
+       compute_error ();
+       output_results (cycle);
+      }
+  }
 }
 
 
@@ -1104,6 +1107,9 @@ int main ()
 
   try
     {
+      using namespace dealii;
+      using namespace Step47;
+
       deallog.depth_console (0);
 
       LaplaceProblem<2> laplace_problem_2d;

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.