]> https://gitweb.dealii.org/ - dealii.git/commitdiff
removed unfinished steps (release tasks step 8). 2157/head
authorTimo Heister <timo.heister@gmail.com>
Thu, 11 Feb 2016 13:08:11 +0000 (08:08 -0500)
committerTimo Heister <timo.heister@gmail.com>
Thu, 11 Feb 2016 13:08:11 +0000 (08:08 -0500)
16 files changed:
examples/step-47/CMakeLists.txt [deleted file]
examples/step-47/doc/builds-on [deleted file]
examples/step-47/doc/intro.dox [deleted file]
examples/step-47/doc/kind [deleted file]
examples/step-47/doc/points.dat [deleted file]
examples/step-47/doc/results.dox [deleted file]
examples/step-47/doc/tooltip [deleted file]
examples/step-47/doc/vertices.dat [deleted file]
examples/step-47/step-47.cc [deleted file]
examples/step-50/CMakeLists.txt [deleted file]
examples/step-50/doc/builds-on [deleted file]
examples/step-50/doc/intro.dox [deleted file]
examples/step-50/doc/kind [deleted file]
examples/step-50/doc/results.dox [deleted file]
examples/step-50/doc/tooltip [deleted file]
examples/step-50/step-50.cc [deleted file]

diff --git a/examples/step-47/CMakeLists.txt b/examples/step-47/CMakeLists.txt
deleted file mode 100644 (file)
index 454844f..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-##
-#  CMake script for the step-47 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "step-47")
-
-# Declare all source files the target consists of. Here, this is only
-# the one step-X.cc file, but as you expand your project you may wish
-# to add other source files as well. If your project becomes much larger,
-# you may want to either replace the following statement by something like
-#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
-#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
-#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
-# or switch altogether to the large project CMakeLists.txt file discussed
-# in the "CMake in user projects" page accessible from the "User info"
-# page of the documentation.
-SET(TARGET_SRC
-  ${TARGET}.cc
-  )
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.4 QUIET
-  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
-  )
-IF(NOT ${deal.II_FOUND})
-  MESSAGE(FATAL_ERROR "\n"
-    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
-    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
-    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
-    )
-ENDIF()
-
-DEAL_II_INITIALIZE_CACHED_VARIABLES()
-PROJECT(${TARGET})
-DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-47/doc/builds-on b/examples/step-47/doc/builds-on
deleted file mode 100644 (file)
index 9fdd726..0000000
+++ /dev/null
@@ -1 +0,0 @@
-step-5
diff --git a/examples/step-47/doc/intro.dox b/examples/step-47/doc/intro.dox
deleted file mode 100644 (file)
index 1518606..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
diff --git a/examples/step-47/doc/kind b/examples/step-47/doc/kind
deleted file mode 100644 (file)
index 6816e90..0000000
+++ /dev/null
@@ -1 +0,0 @@
-unfinished
diff --git a/examples/step-47/doc/points.dat b/examples/step-47/doc/points.dat
deleted file mode 100644 (file)
index 01ce16a..0000000
+++ /dev/null
@@ -1,21 +0,0 @@
-#xfem quadrature Points
-0.781998 0.781998
-0.941587 0.72508
-0.72508 0.941587
-0.926335 0.926335
-0.397329 0.926335
-0.583333 0.72508
-0.583333 0.941587
-0.686004 0.781998
-0.72508 0.583333
-0.926335 0.397329
-0.781998 0.686004
-0.941587 0.583333
-0.115331 0.315492
-0.430422 0.60008
-0.069578 0.816587
-0.259669 0.892842
-0.315492 0.115331
-0.816587 0.069578
-0.60008 0.430422
-0.892842 0.259669
diff --git a/examples/step-47/doc/results.dox b/examples/step-47/doc/results.dox
deleted file mode 100644 (file)
index f4c6fee..0000000
+++ /dev/null
@@ -1 +0,0 @@
-<h1>Results</h1>
diff --git a/examples/step-47/doc/tooltip b/examples/step-47/doc/tooltip
deleted file mode 100644 (file)
index f2003d7..0000000
+++ /dev/null
@@ -1 +0,0 @@
-X-FEM
diff --git a/examples/step-47/doc/vertices.dat b/examples/step-47/doc/vertices.dat
deleted file mode 100644 (file)
index 8d1ef87..0000000
+++ /dev/null
@@ -1,20 +0,0 @@
-#vertices of xfem subcells
-0   0
-1   0
-1   1
-0   1
-
-0.25   1
-1   0.25
-
-0.625   0.625
-0.75   0.75
-
-0.75   0.75
-0.625   1
-
-0.75   0.75
-1   0.625
-
-0   0
-0.625   0.625
diff --git a/examples/step-47/step-47.cc b/examples/step-47/step-47.cc
deleted file mode 100644 (file)
index 6959040..0000000
+++ /dev/null
@@ -1,1137 +0,0 @@
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2011 - 2015 by the deal.II authors
- *
- * This file is part of the deal.II library.
- *
- * The deal.II library is free software; you can use it, redistribute
- * it, and/or modify it under the terms of the GNU Lesser General
- * Public License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- * The full text of the license can be found in the file LICENSE at
- * the top level of the deal.II distribution.
- *
- * ---------------------------------------------------------------------
-
- *
- * Author: Wolfgang Bangerth, University of Heidelberg, 2000
- */
-
-
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/manifold_lib.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_nothing.h>
-#include <deal.II/fe/fe_system.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/numerics/matrix_tools.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <fstream>
-#include <iostream>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/grid/grid_out.h>
-
-
-#include <deal.II/lac/constraint_matrix.h>
-
-#include <deal.II/grid/grid_refinement.h>
-
-#include <deal.II/numerics/error_estimator.h>
-
-namespace Step47
-{
-  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
-  {
-  public:
-    LaplaceProblem ();
-    ~LaplaceProblem ();
-
-    void run ();
-
-  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;
-
-    ConstraintMatrix      constraints;
-
-    SparsityPattern       sparsity_pattern;
-    SparseMatrix<double>  system_matrix;
-
-    Vector<double>        solution;
-    Vector<double>        system_rhs;
-  };
-
-
-
-
-  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;
-
-    virtual void value_list (const std::vector<Point<dim> > &points,
-                             std::vector<double>            &values,
-                             const unsigned int              component = 0) const;
-  };
-
-
-
-  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>
-  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();
-
-    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 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>
-  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>
-  LaplaceProblem<dim>::~LaplaceProblem ()
-  {
-    dof_handler.clear ();
-  }
-
-
-
-  template <int dim>
-  double
-  level_set (const Point<dim> &p)
-  {
-    return p.norm() - 0.5;
-  }
-
-
-
-  template <int dim>
-  Tensor<1,dim>
-  grad_level_set (const Point<dim> &p)
-  {
-    return p / p.norm();
-  }
-
-
-
-  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>
-  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);
-
-    dof_handler.distribute_dofs (fe_collection);
-
-    solution.reinit (dof_handler.n_dofs());
-    system_rhs.reinit (dof_handler.n_dofs());
-
-
-    constraints.clear ();
-//TODO: fix this, it currently crashes
-    // DoFTools::make_hanging_node_constraints (dof_handler, constraints);
-
-//TODO: component 1 must satisfy zero boundary conditions
-    constraints.close();
-
-
-    DynamicSparsityPattern dsp(dof_handler.n_dofs());
-    DoFTools::make_sparsity_pattern (dof_handler, dsp);
-
-    constraints.condense (dsp);
-
-    sparsity_pattern.copy_from(dsp);
-
-    system_matrix.reinit (sparsity_pattern);
-  }
-
-
-  template <int dim>
-  void LaplaceProblem<dim>::assemble_system ()
-  {
-    const QGauss<dim>  quadrature_formula(3);
-
-
-    FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
-                                   update_values    |  update_gradients |
-                                   update_quadrature_points  |  update_JxW_values);
-
-    const unsigned int   n_q_points    = quadrature_formula.size();
-
-    FullMatrix<double>   cell_matrix;
-    Vector<double>       cell_rhs;
-
-    std::vector<types::global_dof_index> local_dof_indices;
-
-    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();
-
-    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;
-
-        if (cell->active_fe_index() == 0)
-          {
-            plain_fe_values.reinit (cell);
-
-            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<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));
-
-
-                  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
-            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);
-      }
-
-
-    std::map<types::global_dof_index,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 the elements on both sides of the discontinuity. The
-// discontinuity line is approximated by a piece-wise linear interpolation
-// between the intersection of the discontinuity with the edges of the
-// elements. The vector level_set_values has the values of the level set
-// function at the vertices of the elements. From these values can be found by
-// linear interpolation the intersections. There are three kind of
-// decomposition that 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;
-      }
-
-    // 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 following three classes of decomposition type 1:
-    // ++++, ---- type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+ type
-    // 3: +--+, ++--, +-+-, -++-, --++, -+-+
-
-    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;
-
-    unsigned int Pos = 100;
-
-    Point<dim> v0(0,0);
-    Point<dim> v1(1,0);
-    Point<dim> v2(0,1);
-    Point<dim> v3(1,1);
-
-    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);
-              }
-          }
-
-        Assert (xfem_quadrature.size() == plain_quadrature.size() * 5, ExcInternalError());
-        return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
-      }
-
-    // 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
-
-        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);
-      }
-
-    return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
-
-  }
-
-  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)
-
-  {
-    // 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.
-
-    unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
-
-    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);
-
-    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;
-          }
-
-          default:
-            Assert (false, ExcNotImplemented());
-          }
-
-
-        Tensor<2,dim> jacobian;
-
-        // 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);
-
-        // 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>
-  void LaplaceProblem<dim>::solve ()
-  {
-    SolverControl           solver_control (1000, 1e-12);
-    SolverCG<>              solver (solver_control);
-
-    PreconditionSSOR<> preconditioner;
-    preconditioner.initialize(system_matrix, 1.2);
-
-    solver.solve (system_matrix, solution, system_rhs,
-                  preconditioner);
-
-    constraints.distribute (solution);
-  }
-
-
-
-  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);
-
-    triangulation.execute_coarsening_and_refinement ();
-  }
-
-
-
-  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;
-
-    virtual std::vector<std::string> get_names () const;
-
-    virtual
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    get_data_component_interpretation () const;
-
-    virtual UpdateFlags get_needed_update_flags () const;
-  };
-
-
-  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;
-  }
-
-
-  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());
-
-    for (unsigned int q=0; q<n_quadrature_points; ++q)
-      {
-        computed_quantities[q](0)
-          = (uh[q](0)
-             +
-//TODO: shift in weight function is missing!
-             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";
-
-    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);
-  }
-
-
-
-  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;
-  }
-
-
-
-
-  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);
-
-            static const SphericalManifold<dim> boundary;
-            triangulation.set_all_manifold_ids_on_boundary(0);
-            triangulation.set_manifold (0, boundary);
-
-            triangulation.refine_global (2);
-          }
-        else
-          triangulation.refine_global (1);
-//      refine_grid ();
-
-
-        std::cout << "   Number of active cells:       "
-                  << triangulation.n_active_cells()
-                  << std::endl;
-
-        setup_system ();
-
-        std::cout << "   Number of degrees of freedom: "
-                  << dof_handler.n_dofs()
-                  << std::endl;
-
-        assemble_system ();
-        solve ();
-        compute_error ();
-        output_results (cycle);
-      }
-  }
-}
-
-
-
-int main ()
-{
-
-  try
-    {
-      using namespace dealii;
-      using namespace Step47;
-
-      LaplaceProblem<2> laplace_problem_2d;
-      laplace_problem_2d.run ();
-    }
-  catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-                << exc.what() << std::endl
-                << "Aborting!" << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-
-      return 1;
-    }
-  catch (...)
-    {
-      std::cerr << std::endl << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-                << "Aborting!" << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      return 1;
-    }
-
-  return 0;
-}
diff --git a/examples/step-50/CMakeLists.txt b/examples/step-50/CMakeLists.txt
deleted file mode 100644 (file)
index 119c115..0000000
+++ /dev/null
@@ -1,57 +0,0 @@
-##
-#  CMake script for the step-32 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "step-50")
-
-# Declare all source files the target consists of. Here, this is only
-# the one step-X.cc file, but as you expand your project you may wish
-# to add other source files as well. If your project becomes much larger,
-# you may want to either replace the following statement by something like
-#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
-#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
-#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
-# or switch altogether to the large project CMakeLists.txt file discussed
-# in the "CMake in user projects" page accessible from the "User info"
-# page of the documentation.
-SET(TARGET_SRC
-  ${TARGET}.cc
-  )
-
-# Define the output that should be cleaned:
-SET(CLEAN_UP_FILES *.vtu *.pvtu *.visit)
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.4 QUIET
-  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
-  )
-IF(NOT ${deal.II_FOUND})
-  MESSAGE(FATAL_ERROR "\n"
-    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
-    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
-    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
-    )
-ENDIF()
-
-#
-# Are all dependencies fulfilled?
-#
-IF( NOT DEAL_II_WITH_MPI OR
-    NOT DEAL_II_WITH_P4EST OR
-    NOT DEAL_II_WITH_TRILINOS )
-  MESSAGE(FATAL_ERROR "
-Error! The deal.II library found at ${DEAL_II_PATH} was not configured with
-    DEAL_II_WITH_MPI = ON
-    DEAL_II_WITH_P4EST = ON
-    DEAL_II_WITH_TRILINOS = ON
-One or all of these are OFF in your installation but are required for this tutorial step."
-    )
-ENDIF()
-
-DEAL_II_INITIALIZE_CACHED_VARIABLES()
-PROJECT(${TARGET})
-DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-50/doc/builds-on b/examples/step-50/doc/builds-on
deleted file mode 100644 (file)
index 79df0ee..0000000
+++ /dev/null
@@ -1 +0,0 @@
-step-16 step-40
diff --git a/examples/step-50/doc/intro.dox b/examples/step-50/doc/intro.dox
deleted file mode 100644 (file)
index 9966672..0000000
+++ /dev/null
@@ -1,90 +0,0 @@
-<br>
-
-<i>This program has evolved from a version originally written by Guido
-Kanschat in 2003. It has undergone significant revisions by B&auml;rbel
-Janssen, Guido Kanschat and Wolfgang Bangerth in 2009 and 2010 to demonstrate
-multigrid algorithms on adaptively refined meshes.
-</i>
-
-
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-
-This example shows the basic usage of the multilevel functions in
-deal.II. It solves the same problem as used in step-6,
-but demonstrating the things one has to provide when using multigrid
-as a preconditioner. In particular, this requires that we define a
-hierarchy of levels, provide transfer operators from one level to the
-next and back, and provide representations of the Laplace operator on
-each level.
-
-In order to allow sufficient flexibility in conjunction with systems of
-differential equations and block preconditioners, quite a few different objects
-have to be created before starting the multilevel method, although
-most of what needs to be done is provided by deal.II itself. These are
-<ul>
-<li>An the object handling transfer between grids; we use the
-    MGTransferPrebuilt class for this that does almost all of the work
-    inside the library.
-<li>The solver on the coarsest level; here, we use MGCoarseGridHouseholder.
-<li>The smoother on all other levels, which in our case will be the
-    MGSmootherRelaxation class using SOR as the underlying method
-<li>And MGMatrix, a class having a special level multiplication, i.e. we
-    basically store one matrix per grid level and allow multiplication
-    with it.
-</ul>
-Most of these objects will only be needed inside the function that
-actually solves the linear system. There, these objects are combined
-in an object of type Multigrid, containing the implementation of the
-V-cycle, which is in turn used by the preconditioner PreconditionMG,
-ready for plug-in into a linear solver of the LAC library.
-
-The multilevel method in deal.II follows in many respects the outlines
-of the various publications by James Bramble, Joseph Pasciak and
-Jinchao Xu (i.e. the "BPX" framework). In order to understand many of
-the options, a rough familiarity with their work is quite helpful.
-
-However, in comparison to this framework, the implementation in
-deal.II has to take into account the fact that we want to solve linear
-systems on adaptively refined meshes. This leads to the complication
-that it isn't quite as clear any more what exactly a "level" in a
-multilevel hierarchy of a mesh is. The following image shows what we
-consider to be a "level":
-
-<p align="center">
-  @image html "hanging_nodes.png" ""
-</p>
-
-In other words, the fine level in this mesh consists only of the
-degrees of freedom that are defined on the refined cells, but does not
-extend to that part of the domain that is not refined. While this
-guarantees that the overall effort grows as ${\cal O}(N)$ as necessary
-for optimal multigrid complexity, it leads to problems when defining
-where to smooth and what boundary conditions to pose for the operators
-defined on individual levels if the level boundary is not an external
-boundary. These questions are discussed in detail in the
-@ref mg_paper "Multigrid paper by Janssen and Kanschat" that describes
-the implementation in deal.II.
-
-
-
-<h3>The testcase</h3>
-
-The problem we solve here is exactly the same as in
-step-6, the only difference being the solver we use
-here. You may want to look there for a definition of what we solve,
-right hand side and boundary conditions. Obviously, the program would
-also work if we changed the geometry and other pieces of data that
-defines this particular problem.
-
-The things that are new are all those parts that concern the
-multigrid. In particular, this includes the following members of the
-main class:
-- <code>LaplaceProblem::mg_dof_handler</code>
-- <code>LaplaceProblem::mg_sparsity</code>
-- <code>LaplaceProblem::mg_matrices</code>
-- <code>LaplaceProblem::mg_interface_matrices_up</code>
-- <code>LaplaceProblem::assemble_multigrid ()</code>
-- <code>LaplaceProblem::solve ()</code>
-Take a look at these functions.
diff --git a/examples/step-50/doc/kind b/examples/step-50/doc/kind
deleted file mode 100644 (file)
index 6816e90..0000000
+++ /dev/null
@@ -1 +0,0 @@
-unfinished
diff --git a/examples/step-50/doc/results.dox b/examples/step-50/doc/results.dox
deleted file mode 100644 (file)
index 2d5fe0d..0000000
+++ /dev/null
@@ -1,98 +0,0 @@
-<h1>Results</h1>
-
-The output that this program generates is, of course, the same as that
-of step-6, so you may see there for more results. On the
-other hand, since no tutorial program is a good one unless it has at
-least one colorful picture, here is, again, the solution:
-
-
-When run, the output of this program is
-<pre>
-Cycle 0:
-   Number of active cells:       20
-   Number of degrees of freedom: 25 (by level: 8, 25)
-   7 CG iterations needed to obtain convergence.
-Cycle 1:
-   Number of active cells:       44
-   Number of degrees of freedom: 57 (by level: 8, 25, 48)
-   8 CG iterations needed to obtain convergence.
-Cycle 2:
-   Number of active cells:       92
-   Number of degrees of freedom: 117 (by level: 8, 25, 80, 60)
-   9 CG iterations needed to obtain convergence.
-Cycle 3:
-   Number of active cells:       188
-   Number of degrees of freedom: 221 (by level: 8, 25, 80, 200)
-   12 CG iterations needed to obtain convergence.
-Cycle 4:
-   Number of active cells:       416
-   Number of degrees of freedom: 485 (by level: 8, 25, 89, 288, 280)
-   13 CG iterations needed to obtain convergence.
-Cycle 5:
-   Number of active cells:       800
-   Number of degrees of freedom: 925 (by level: 8, 25, 89, 288, 784, 132)
-   14 CG iterations needed to obtain convergence.
-Cycle 6:
-   Number of active cells:       1628
-   Number of degrees of freedom: 1865 (by level: 8, 25, 89, 304, 1000, 1164, 72)
-   14 CG iterations needed to obtain convergence.
-Cycle 7:
-   Number of active cells:       3194
-   Number of degrees of freedom: 3603 (by level: 8, 25, 89, 328, 1032, 2200, 1392)
-   16 CG iterations needed to obtain convergence.
-</pre>
-That's not perfect &mdash; we would have hoped for a constant number
-of iterations rather than one that increases as we get more and more
-degrees of freedom &mdash; but it is also not far away. The reason for
-this is easy enough to understand, however: since we have a strongly
-varying coefficient, the operators that we assembly by quadrature on
-the lower levels become worse and worse approximations of the operator
-on the finest level. Consequently, even if we had perfect solvers on
-the coarser levels, they would not be good preconditioners on the
-finest level. This theory is easily tested by comparing results when
-we use a constant coefficient: in that case, the number of iterations
-remains constant at 9 after the first three or four refinement steps.
-
-We can also compare what this program produces with how @ref step_5
-"step-5" performed. To solve the same problem as in step-5, the only
-two changes that are necessary are (i) to replace the body of the
-function <code>LaplaceProblem::refine_grid</code> by a call to
-<code>triangulation.refine_global(1)</code>, and (ii) to use the same
-SolverControl object and tolerance as in step-5 &mdash; the rest of the
-program remains unchanged. In that case, here is how the solvers used
-in step-5 and the multigrid solver used in the current program
-compare:
-<table align="center">
-<tr><th>cells</th><th>step-5</th><th>step-16</th></tr>
-<tr><td>20</td>   <td>13</td> <td>6</td> </tr>
-<tr><td>80</td>   <td>17</td> <td>7</td> </tr>
-<tr><td>320</td>  <td>29</td> <td>9</td> </tr>
-<tr><td>1280</td> <td>51</td> <td>10</td> </tr>
-<tr><td>5120</td> <td>94</td> <td>11</td> </tr>
-<tr><td>20480</td><td>180</td><td>13</td></tr>
-</table>
-This isn't only fewer iterations than in step-5 (each of which
-is, however, much more expensive) but more importantly, the number of
-iterations also grows much more slowly under mesh refinement (again,
-it would be almost constant if the coefficient was constant rather
-than strongly varying as chosen here). This justifies the common
-observation that, whenever possible, multigrid methods should be used
-for second order problems.
-
-
-<h3> Possible extensions </h3>
-
-A close inspection of this program's performance shows that it is mostly
-dominated by matrix-vector operations. step-37 shows one way
-how this can be avoided by working with matrix-free methods.
-
-Another avenue would be to use algebraic multigrid methods. The
-geometric multigrid method used here can at times be a bit awkward to
-implement because it needs all those additional data structures, and
-it becomes even more difficult if the program is to run in %parallel on
-machines coupled through MPI, for example. In that case, it would be
-simpler if one could use a black-box preconditioner that uses some
-sort of multigrid hierarchy for good performance but can figure out
-level matrices and similar things out by itself. Algebraic multigrid
-methods do exactly this, and we will use them in
-step-31 for the solution of a Stokes problem.
diff --git a/examples/step-50/doc/tooltip b/examples/step-50/doc/tooltip
deleted file mode 100644 (file)
index 2fd6559..0000000
+++ /dev/null
@@ -1 +0,0 @@
-Multigrid on adaptive meshes.
diff --git a/examples/step-50/step-50.cc b/examples/step-50/step-50.cc
deleted file mode 100644 (file)
index a914980..0000000
+++ /dev/null
@@ -1,1021 +0,0 @@
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2003 - 2015 by the deal.II authors
- *
- * This file is part of the deal.II library.
- *
- * The deal.II library is free software; you can use it, redistribute
- * it, and/or modify it under the terms of the GNU Lesser General
- * Public License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- * The full text of the license can be found in the file LICENSE at
- * the top level of the deal.II distribution.
- *
- * ---------------------------------------------------------------------
-
- *
- * Author: Guido Kanschat and Timo Heister
- */
-
-
-// @note: This a work in progress example of parallel geometric
-// multigrid. Some parts are still in heavy development.
-
-// This program is a parallel version of step-16 with a slightly different
-// problem setup.
-
-// @sect3{Include files}
-
-// Again, the first few include files
-// are already known, so we won't
-// comment on them:
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/base/conditional_ostream.h>
-
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/lac/precondition.h>
-
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-
-#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/error_estimator.h>
-
-#include <deal.II/base/index_set.h>
-#include <deal.II/distributed/tria.h>
-#include <deal.II/distributed/grid_refinement.h>
-
-#include <deal.II/multigrid/mg_constrained_dofs.h>
-#include <deal.II/multigrid/multigrid.h>
-#include <deal.II/multigrid/mg_transfer.h>
-#include <deal.II/multigrid/mg_tools.h>
-#include <deal.II/multigrid/mg_coarse.h>
-#include <deal.II/multigrid/mg_smoother.h>
-#include <deal.II/multigrid/mg_matrix.h>
-
-
-#include <deal.II/lac/generic_linear_algebra.h>
-
-// #define USE_PETSC_LA PETSc is not quite supported yet
-
-namespace LA
-{
-#ifdef USE_PETSC_LA
-  using namespace dealii::LinearAlgebraPETSc;
-#else
-  using namespace dealii::LinearAlgebraTrilinos;
-#endif
-}
-
-// This is C++:
-#include <iostream>
-#include <fstream>
-#include <sstream>
-
-// The last step is as in all
-// previous programs:
-namespace Step50
-{
-  using namespace dealii;
-
-
-  // @sect3{The <code>LaplaceProblem</code> class template}
-
-  // This main class is very similar to step-16, except that we are storing a
-  // parallel Triangulation and parallel versions of matrices and vectors.
-  template <int dim>
-  class LaplaceProblem
-  {
-  public:
-    LaplaceProblem (const unsigned int deg);
-    void run ();
-
-  private:
-    void setup_system ();
-    void assemble_system ();
-    void assemble_multigrid ();
-    void solve ();
-    void refine_grid ();
-    void output_results (const unsigned int cycle) const;
-
-    ConditionalOStream                        pcout;
-
-    parallel::distributed::Triangulation<dim>   triangulation;
-    FE_Q<dim>            fe;
-    DoFHandler<dim>    mg_dof_handler;
-
-    typedef LA::MPI::SparseMatrix matrix_t;
-    typedef LA::MPI::Vector vector_t;
-
-    matrix_t system_matrix;
-
-    IndexSet locally_relevant_set;
-
-    ConstraintMatrix     hanging_node_constraints;
-    ConstraintMatrix     constraints;
-
-    vector_t       solution;
-    vector_t       system_rhs;
-
-    const unsigned int degree;
-
-    // Finally we are storing the various parallel multigrid matrices. Our
-    // problem is self-adjoint, so the interface matrices are the transpose
-    // of each other, so we only need to compute/store them once.
-    MGLevelObject<matrix_t> mg_matrices;
-    MGLevelObject<matrix_t> mg_interface_matrices;
-    //
-    MGConstrainedDoFs                    mg_constrained_dofs;
-  };
-
-
-
-  // @sect3{Nonconstant coefficients}
-
-  // The implementation of nonconstant
-  // coefficients is copied verbatim
-  // from step-5 and step-6:
-
-  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;
-
-    virtual void value_list (const std::vector<Point<dim> > &points,
-                             std::vector<double>            &values,
-                             const unsigned int              component = 0) const;
-  };
-
-
-
-  template <int dim>
-  double Coefficient<dim>::value (const Point<dim> &p,
-                                  const unsigned int) const
-  {
-    if (p.square() < 0.5*0.5)
-      return 5;
-    else
-      return 1;
-  }
-
-
-
-  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();
-
-    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)
-      values[i] = Coefficient<dim>::value (points[i]);
-  }
-
-
-  // @sect3{The <code>LaplaceProblem</code> class implementation}
-
-  // @sect4{LaplaceProblem::LaplaceProblem}
-
-  // The constructor is left mostly
-  // unchanged. We take the polynomial degree
-  // of the finite elements to be used as a
-  // constructor argument and store it in a
-  // member variable.
-  //
-  // By convention, all adaptively refined
-  // triangulations in deal.II never change by
-  // more than one level across a face between
-  // cells. For our multigrid algorithms,
-  // however, we need a slightly stricter
-  // guarantee, namely that the mesh also does
-  // not change by more than refinement level
-  // across vertices that might connect two
-  // cells. In other words, we must prevent the
-  // following situation:
-  //
-  // @image html limit_level_difference_at_vertices.png ""
-  //
-  // This is achieved by passing the
-  // Triangulation::limit_level_difference_at_vertices
-  // flag to the constructor of the
-  // triangulation class.
-  template <int dim>
-  LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
-    :
-    pcout (std::cout,
-           (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)
-            == 0)),
-    triangulation (MPI_COMM_WORLD,Triangulation<dim>::
-                   limit_level_difference_at_vertices,
-                   parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
-    fe (degree),
-    mg_dof_handler (triangulation),
-    degree(degree)
-  {}
-
-
-  // @sect4{LaplaceProblem::setup_system}
-
-  // The following function extends what the
-  // corresponding one in step-6 did. The top
-  // part, apart from the additional output,
-  // does the same:
-  template <int dim>
-  void LaplaceProblem<dim>::setup_system ()
-  {
-    mg_dof_handler.distribute_dofs (fe);
-    mg_dof_handler.distribute_mg_dofs (fe);
-
-    DoFTools::extract_locally_relevant_dofs (mg_dof_handler,
-                                             locally_relevant_set);
-
-    solution.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
-    system_rhs.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
-
-    // But it starts to be a wee bit different
-    // here, although this still doesn't have
-    // anything to do with multigrid
-    // methods. step-6 took care of boundary
-    // values and hanging nodes in a separate
-    // step after assembling the global matrix
-    // from local contributions. This works,
-    // but the same can be done in a slightly
-    // simpler way if we already take care of
-    // these constraints at the time of copying
-    // local contributions into the global
-    // matrix. To this end, we here do not just
-    // compute the constraints do to hanging
-    // nodes, but also due to zero boundary
-    // conditions. We will
-    // use this set of constraints later on to
-    // help us copy local contributions
-    // correctly into the global linear system
-    // right away, without the need for a later
-    // clean-up stage:
-    constraints.reinit (locally_relevant_set);
-    hanging_node_constraints.reinit (locally_relevant_set);
-    DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints);
-    DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
-
-    typename FunctionMap<dim>::type      dirichlet_boundary;
-    ConstantFunction<dim>                    homogeneous_dirichlet_bc (1.0);
-    dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
-    VectorTools::interpolate_boundary_values (mg_dof_handler,
-                                              dirichlet_boundary,
-                                              constraints);
-    constraints.close ();
-    hanging_node_constraints.close ();
-
-    DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(), mg_dof_handler.n_dofs());
-    DoFTools::make_sparsity_pattern (mg_dof_handler, dsp, constraints);
-    system_matrix.reinit (mg_dof_handler.locally_owned_dofs(), dsp, MPI_COMM_WORLD, true);
-
-
-    // The multigrid constraints have to be
-    // initialized. They need to know about
-    // the boundary values as well, so we
-    // pass the <code>dirichlet_boundary</code>
-    // here as well.
-    mg_constrained_dofs.clear();
-    mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary);
-
-
-    // Now for the things that concern the
-    // multigrid data structures. First, we
-    // resize the multilevel objects to hold
-    // matrices and sparsity patterns for every
-    // level. The coarse level is zero (this is
-    // mandatory right now but may change in a
-    // future revision). Note that these
-    // functions take a complete, inclusive
-    // range here (not a starting index and
-    // size), so the finest level is
-    // <code>n_levels-1</code>.  We first have
-    // to resize the container holding the
-    // SparseMatrix classes, since they have to
-    // release their SparsityPattern before the
-    // can be destroyed upon resizing.
-    const unsigned int n_levels = triangulation.n_global_levels();
-
-    mg_interface_matrices.resize(0, n_levels-1);
-    mg_interface_matrices.clear ();
-    mg_matrices.resize(0, n_levels-1);
-    mg_matrices.clear ();
-
-    // Now, we have to provide a matrix on each
-    // level. To this end, we first use the
-    // MGTools::make_sparsity_pattern function
-    // to first generate a preliminary
-    // compressed sparsity pattern on each
-    // level (see the @ref Sparsity module for
-    // more information on this topic) and then
-    // copy it over to the one we really
-    // want. The next step is to initialize
-    // both kinds of level matrices with these
-    // sparsity patterns.
-    //
-    // It may be worth pointing out that the
-    // interface matrices only have entries for
-    // degrees of freedom that sit at or next
-    // to the interface between coarser and
-    // finer levels of the mesh. They are
-    // therefore even sparser than the matrices
-    // on the individual levels of our
-    // multigrid hierarchy. If we were more
-    // concerned about memory usage (and
-    // possibly the speed with which we can
-    // multiply with these matrices), we should
-    // use separate and different sparsity
-    // patterns for these two kinds of
-    // matrices.
-    for (unsigned int level=0; level<n_levels; ++level)
-      {
-        DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(level),
-                                   mg_dof_handler.n_dofs(level));
-        MGTools::make_sparsity_pattern(mg_dof_handler, dsp, level);
-
-        mg_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
-                                  mg_dof_handler.locally_owned_mg_dofs(level),
-                                  dsp,
-                                  MPI_COMM_WORLD, true);
-
-        mg_interface_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
-                                            mg_dof_handler.locally_owned_mg_dofs(level),
-                                            dsp,
-                                            MPI_COMM_WORLD, true);
-      }
-  }
-
-
-  // @sect4{LaplaceProblem::assemble_system}
-
-  // The following function assembles the
-  // linear system on the finest level of the
-  // mesh. It is almost exactly the same as in
-  // step-6, with the exception that we don't
-  // eliminate hanging nodes and boundary
-  // values after assembling, but while copying
-  // local contributions into the global
-  // matrix. This is not only simpler but also
-  // more efficient for large problems.
-  //
-  // This latter trick is something that only
-  // found its way into deal.II over time and
-  // wasn't used in the initial version of this
-  // tutorial program. There is, however, a
-  // discussion of this function in the
-  // introduction of step-27.
-  template <int dim>
-  void LaplaceProblem<dim>::assemble_system ()
-  {
-    const QGauss<dim>  quadrature_formula(degree+1);
-
-    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<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    const Coefficient<dim> coefficient;
-    std::vector<double>    coefficient_values (n_q_points);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = mg_dof_handler.begin_active(),
-    endc = mg_dof_handler.end();
-    for (; cell!=endc; ++cell)
-      if (cell->is_locally_owned())
-        {
-          cell_matrix = 0;
-          cell_rhs = 0;
-
-          fe_values.reinit (cell);
-
-          coefficient.value_list (fe_values.get_quadrature_points(),
-                                  coefficient_values);
-
-          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] *
-                                       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) *
-                                10.0 *
-                                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);
-        }
-
-    system_matrix.compress(VectorOperation::add);
-    system_rhs.compress(VectorOperation::add);
-  }
-
-
-  // @sect4{LaplaceProblem::assemble_multigrid}
-
-  // The next function is the one that builds
-  // the linear operators (matrices) that
-  // define the multigrid method on each level
-  // of the mesh. The integration core is the
-  // same as above, but the loop below will go
-  // over all existing cells instead of just
-  // the active ones, and the results must be
-  // entered into the correct matrix. Note also
-  // that since we only do multilevel
-  // preconditioning, no right-hand side needs
-  // to be assembled here.
-  //
-  // Before we go there, however, we have to
-  // take care of a significant amount of book
-  // keeping:
-  template <int dim>
-  void LaplaceProblem<dim>::assemble_multigrid ()
-  {
-    QGauss<dim>  quadrature_formula(1+degree);
-
-    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);
-
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    const Coefficient<dim> coefficient;
-    std::vector<double>    coefficient_values (n_q_points);
-
-    // Next a few things that are specific to building the multigrid
-    // data structures (since we only need them in the current
-    // function, rather than also elsewhere, we build them here
-    // instead of the <code>setup_system</code> function). Some of the
-    // following may be a bit obscure if you're not familiar with the
-    // algorithm actually implemented in deal.II to support multilevel
-    // algorithms on adaptive meshes; if some of the things below seem
-    // strange, take a look at the @ref mg_paper.
-    //
-    // Our first job is to identify those degrees of freedom on each level
-    // that are located on interfaces between adaptively refined levels, and
-    // those that lie on the interface but also on the exterior boundary of
-    // the domain. The <code>MGConstraints</code> already computed the
-    // information for us when we called initialize in
-
-    // <code>setup_system()</code>.
-    // of type IndexSet on each level (get_refinement_edge_indices(),
-
-    // The indices just identified will later be used to decide where
-    // the assembled value has to be added into on each level.  On the
-    // other hand, we also have to impose zero boundary conditions on
-    // the external boundary of each level. But this the
-    // <code>MGConstraints</code> knows it. So we simply ask for them
-    // by calling <code>get_boundary_indices ()</code>.  The third
-    // step is to construct constraints on all those degrees of
-    // freedom: their value should be zero after each application of
-    // the level operators. To this end, we construct ConstraintMatrix
-    // objects for each level, and add to each of these constraints
-    // for each degree of freedom. Due to the way the ConstraintMatrix
-    // stores its data, the function to add a constraint on a single
-    // degree of freedom and force it to be zero is called
-    // ConstraintMatrix::add_line(); doing so for several degrees of
-    // freedom at once can be done using
-    // ConstraintMatrix::add_lines():
-    std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_global_levels());
-    ConstraintMatrix empty_constraints;
-    for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
-      {
-        IndexSet dofset;
-        DoFTools::extract_locally_relevant_level_dofs (mg_dof_handler, level, dofset);
-        boundary_constraints[level].reinit(dofset);
-        boundary_constraints[level].add_lines (mg_constrained_dofs.get_refinement_edge_indices(level));
-        boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices(level));
-
-        boundary_constraints[level].close ();
-      }
-
-    // Now that we're done with most of our preliminaries, let's start
-    // the integration loop. It looks mostly like the loop in
-    // <code>assemble_system</code>, with two exceptions: (i) we don't
-    // need a right hand side, and more significantly (ii) we don't
-    // just loop over all active cells, but in fact all cells, active
-    // or not. Consequently, the correct iterator to use is
-    // MGDoFHandler::cell_iterator rather than
-    // MGDoFHandler::active_cell_iterator. Let's go about it:
-    typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
-                                            endc = mg_dof_handler.end();
-
-    for (; cell!=endc; ++cell)
-      if (cell->level_subdomain_id()==triangulation.locally_owned_subdomain())
-        {
-          cell_matrix = 0;
-          fe_values.reinit (cell);
-
-          coefficient.value_list (fe_values.get_quadrature_points(),
-                                  coefficient_values);
-
-          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] *
-                                     fe_values.shape_grad(i,q_point) *
-                                     fe_values.shape_grad(j,q_point) *
-                                     fe_values.JxW(q_point));
-
-          // The rest of the assembly is again slightly
-          // different. This starts with a gotcha that is easily
-          // forgotten: The indices of global degrees of freedom we
-          // want here are the ones for current level, not for the
-          // global matrix. We therefore need the function
-          // MGDoFAccessorLLget_mg_dof_indices, not
-          // MGDoFAccessor::get_dof_indices as used in the assembly of
-          // the global system:
-          cell->get_mg_dof_indices (local_dof_indices);
-
-          // Next, we need to copy local contributions into the level
-          // objects. We can do this in the same way as in the global
-          // assembly, using a constraint object that takes care of
-          // constrained degrees (which here are only boundary nodes,
-          // as the individual levels have no hanging node
-          // constraints). Note that the
-          // <code>boundary_constraints</code> object makes sure that
-          // the level matrices contains no contributions from degrees
-          // of freedom at the interface between cells of different
-          // refinement level.
-          boundary_constraints[cell->level()]
-          .distribute_local_to_global (cell_matrix,
-                                       local_dof_indices,
-                                       mg_matrices[cell->level()]);
-
-          // The next step is again slightly more obscure (but
-          // explained in the @ref mg_paper): We need the remainder of
-          // the operator that we just copied into the
-          // <code>mg_matrices</code> object, namely the part on the
-          // interface between cells at the current level and cells
-          // one level coarser. This matrix exists in two directions:
-          // for interior DoFs (index $i$) of the current level to
-          // those sitting on the interface (index $j$), and the other
-          // way around. Of course, since we have a symmetric
-          // operator, one of these matrices is the transpose of the
-          // other.
-          //
-          // The way we assemble these matrices is as follows: since
-          // the are formed from parts of the local contributions, we
-          // first delete all those parts of the local contributions
-          // that we are not interested in, namely all those elements
-          // of the local matrix for which not $i$ is an interface DoF
-          // and $j$ is not. The result is one of the two matrices
-          // that we are interested in, and we then copy it into the
-          // <code>mg_interface_matrices</code> object. The
-          // <code>boundary_interface_constraints</code> object at the
-          // same time makes sure that we delete contributions from
-          // all degrees of freedom that are not only on the interface
-          // but also on the external boundary of the domain.
-          //
-          // The last part to remember is how to get the other
-          // matrix. Since it is only the transpose, we will later (in
-          // the <code>solve()</code> function) be able to just pass
-          // the transpose matrix where necessary.
-
-          const IndexSet &interface_dofs_on_level
-            = mg_constrained_dofs.get_refinement_edge_indices(cell->level());
-          const unsigned int lvl = cell->level();
-
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            for (unsigned int j=0; j<dofs_per_cell; ++j)
-              if (interface_dofs_on_level.is_element(local_dof_indices[i])   // at_refinement_edge(i)
-                  &&
-                  !interface_dofs_on_level.is_element(local_dof_indices[j])   // !at_refinement_edge(j)
-                  &&
-                  (
-                    (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
-                     &&
-                     !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
-                    ) // ( !boundary(i) && !boundary(j) )
-                    ||
-                    (
-                      mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
-                      &&
-                      local_dof_indices[i]==local_dof_indices[j]
-                    ) // ( boundary(i) && boundary(j) && i==j )
-                  )
-                 )
-                {
-                  // do nothing, so add entries to interface matrix
-                }
-              else
-                {
-                  cell_matrix(i,j) = 0;
-                }
-
-
-          empty_constraints
-          .distribute_local_to_global (cell_matrix,
-                                       local_dof_indices,
-                                       mg_interface_matrices[cell->level()]);
-        }
-
-    for (unsigned int i=0; i<triangulation.n_global_levels(); ++i)
-      {
-        mg_matrices[i].compress(VectorOperation::add);
-        mg_interface_matrices[i].compress(VectorOperation::add);
-      }
-  }
-
-
-
-  // @sect4{LaplaceProblem::solve}
-
-  // This is the other function that is significantly different in
-  // support of the multigrid solver (or, in fact, the preconditioner
-  // for which we use the multigrid method).
-  //
-  // Let us start out by setting up two of the components of
-  // multilevel methods: transfer operators between levels, and a
-  // solver on the coarsest level. In finite element methods, the
-  // transfer operators are derived from the finite element function
-  // spaces involved and can often be computed in a generic way
-  // independent of the problem under consideration. In that case, we
-  // can use the MGTransferPrebuilt class that, given the constraints
-  // on the global level and an MGDoFHandler object computes the
-  // matrices corresponding to these transfer operators.
-  //
-  // The second part of the following lines deals with the coarse grid
-  // solver. Since our coarse grid is very coarse indeed, we decide
-  // for a direct solver (a Householder decomposition of the coarsest
-  // level matrix), even if its implementation is not particularly
-  // sophisticated. If our coarse mesh had many more cells than the
-  // five we have here, something better suited would obviously be
-  // necessary here.
-  template <int dim>
-  void LaplaceProblem<dim>::solve ()
-  {
-    // Create the object that deals with the transfer between
-    // different refinement levels. We need to pass it the hanging
-    // node constraints.
-    MGTransferPrebuilt<vector_t> mg_transfer(hanging_node_constraints, mg_constrained_dofs);
-    // Now the prolongation matrix has to be built.  This matrix needs
-    // to take the boundary values on each level into account and
-    // needs to know about the indices at the refinement edges. The
-    // <code>MGConstraints</code> knows about that so pass it as an
-    // argument.
-    mg_transfer.build_matrices(mg_dof_handler);
-
-    matrix_t &coarse_matrix = mg_matrices[0];
-
-    SolverControl coarse_solver_control (1000, 1e-10, false, false);
-    SolverCG<vector_t> coarse_solver(coarse_solver_control);
-    PreconditionIdentity id;
-    MGCoarseGridLACIteration<SolverCG<vector_t>,vector_t> coarse_grid_solver(coarse_solver,
-        coarse_matrix,
-        id);
-
-    // The next component of a multilevel solver or preconditioner is
-    // that we need a smoother on each level. A common choice for this
-    // is to use the application of a relaxation method (such as the
-    // SOR, Jacobi or Richardson method). The MGSmootherPrecondition
-    // class provides support for this kind of smoother. Here, we opt
-    // for the application of a single SOR iteration. To this end, we
-    // define an appropriate <code>typedef</code> and then setup a
-    // smoother object.
-    //
-    // The last step is to initialize the smoother object with our
-    // level matrices and to set some smoothing parameters.  The
-    // <code>initialize()</code> function can optionally take
-    // additional arguments that will be passed to the smoother object
-    // on each level. In the current case for the SOR smoother, this
-    // could, for example, include a relaxation parameter. However, we
-    // here leave these at their default values. The call to
-    // <code>set_steps()</code> indicates that we will use two pre-
-    // and two post-smoothing steps on each level; to use a variable
-    // number of smoother steps on different levels, more options can
-    // be set in the constructor call to the <code>mg_smoother</code>
-    // object.
-    //
-    // The last step results from the fact that
-    // we use the SOR method as a smoother -
-    // which is not symmetric - but we use the
-    // conjugate gradient iteration (which
-    // requires a symmetric preconditioner)
-    // below, we need to let the multilevel
-    // preconditioner make sure that we get a
-    // symmetric operator even for nonsymmetric
-    // smoothers:
-    typedef LA::MPI::PreconditionJacobi Smoother;
-    MGSmootherPrecondition<matrix_t, Smoother, vector_t> mg_smoother;
-    mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5));
-    mg_smoother.set_steps(2);
-    //mg_smoother.set_symmetric(false);
-
-    // The next preparatory step is that we
-    // must wrap our level and interface
-    // matrices in an object having the
-    // required multiplication functions. We
-    // will create two objects for the
-    // interface objects going from coarse to
-    // fine and the other way around; the
-    // multigrid algorithm will later use the
-    // transpose operator for the latter
-    // operation, allowing us to initialize
-    // both up and down versions of the
-    // operator with the matrices we already
-    // built:
-    mg::Matrix<vector_t> mg_matrix(mg_matrices);
-    mg::Matrix<vector_t> mg_interface_up(mg_interface_matrices);
-    mg::Matrix<vector_t> mg_interface_down(mg_interface_matrices);
-
-    // Now, we are ready to set up the
-    // V-cycle operator and the
-    // multilevel preconditioner.
-    Multigrid<vector_t > mg(mg_dof_handler,
-                            mg_matrix,
-                            coarse_grid_solver,
-                            mg_transfer,
-                            mg_smoother,
-                            mg_smoother);
-    //mg.set_debug(6);
-    mg.set_edge_matrices(mg_interface_down, mg_interface_up);
-
-    PreconditionMG<dim, vector_t, MGTransferPrebuilt<vector_t> >
-    preconditioner(mg_dof_handler, mg, mg_transfer);
-
-
-    // With all this together, we can finally
-    // get about solving the linear system in
-    // the usual way:
-    SolverControl solver_control (500, 1e-8*system_rhs.l2_norm(), false);
-    SolverCG<vector_t> solver (solver_control);
-
-    if (false)
-      {
-        /*
-         // code to optionally compare to Trilinos ML
-         TrilinosWrappers::PreconditionAMG prec;
-
-         TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data;
-         //    Amg_data.constant_modes = constant_modes;
-         Amg_data.elliptic = true;
-         Amg_data.higher_order_elements = true;
-         Amg_data.smoother_sweeps = 2;
-         Amg_data.aggregation_threshold = 0.02;
-         // Amg_data.symmetric = true;
-
-         prec.initialize (system_matrix,
-                          Amg_data);
-         solver.solve (system_matrix, solution, system_rhs, prec);
-        */
-      }
-    else
-      {
-        solver.solve (system_matrix, solution, system_rhs,
-                      preconditioner);
-      }
-    pcout << "   CG converged in " << solver_control.last_step() << " iterations." << std::endl;
-
-    constraints.distribute (solution);
-  }
-
-
-
-  // @sect4{Postprocessing}
-
-  // The following two functions postprocess a
-  // solution once it is computed. In
-  // particular, the first one refines the mesh
-  // at the beginning of each cycle while the
-  // second one outputs results at the end of
-  // each such cycle. The functions are almost
-  // unchanged from those in step-6, with the
-  // exception of two minor differences: The
-  // KellyErrorEstimator::estimate function
-  // wants an argument of type DoFHandler, not
-  // MGDoFHandler, and so we have to cast from
-  // derived to base class; and we generate
-  // output in VTK format, to use the more
-  // modern visualization programs available
-  // today compared to those that were
-  // available when step-6 was written.
-  template <int dim>
-  void LaplaceProblem<dim>::refine_grid ()
-  {
-    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-    LA::MPI::Vector temp_solution;
-    temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD);
-    temp_solution = solution;
-
-    KellyErrorEstimator<dim>::estimate (static_cast<DoFHandler<dim>&>(mg_dof_handler),
-                                        QGauss<dim-1>(degree+1),
-                                        typename FunctionMap<dim>::type(),
-                                        temp_solution,
-                                        estimated_error_per_cell);
-
-    parallel::distributed::GridRefinement::
-    refine_and_coarsen_fixed_fraction (triangulation,
-                                       estimated_error_per_cell,
-                                       0.3, 0.0);
-
-    triangulation.execute_coarsening_and_refinement ();
-  }
-
-
-
-  template <int dim>
-  void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
-  {
-    DataOut<dim> data_out;
-
-    LA::MPI::Vector temp_solution;
-    temp_solution.reinit(locally_relevant_set, MPI_COMM_WORLD);
-    temp_solution = solution;
-
-
-    LA::MPI::Vector temp = solution;
-    system_matrix.residual(temp,solution,system_rhs);
-    LA::MPI::Vector res_ghosted = temp_solution;
-    res_ghosted = temp;
-
-    data_out.attach_dof_handler (mg_dof_handler);
-    data_out.add_data_vector (temp_solution, "solution");
-    data_out.add_data_vector (res_ghosted, "res");
-    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 (0);
-
-    const std::string filename = ("solution-" +
-                                  Utilities::int_to_string (cycle, 5) +
-                                  "." +
-                                  Utilities::int_to_string
-                                  (triangulation.locally_owned_subdomain(), 4) +
-                                  ".vtu");
-    std::ofstream output (filename.c_str());
-    data_out.write_vtu (output);
-
-    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
-      {
-        std::vector<std::string> filenames;
-        for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
-          filenames.push_back (std::string("solution-") +
-                               Utilities::int_to_string (cycle, 5) +
-                               "." +
-                               Utilities::int_to_string(i, 4) +
-                               ".vtu");
-        const std::string
-        pvtu_master_filename = ("solution-" +
-                                Utilities::int_to_string (cycle, 5) +
-                                ".pvtu");
-        std::ofstream pvtu_master (pvtu_master_filename.c_str());
-        data_out.write_pvtu_record (pvtu_master, filenames);
-
-        const std::string
-        visit_master_filename = ("solution-" +
-                                 Utilities::int_to_string (cycle, 5) +
-                                 ".visit");
-        std::ofstream visit_master (visit_master_filename.c_str());
-        data_out.write_visit_record (visit_master, filenames);
-
-        std::cout << "   wrote " << pvtu_master_filename << std::endl;
-
-      }
-  }
-
-
-  // @sect4{LaplaceProblem::run}
-
-  // Like several of the functions above, this
-  // is almost exactly a copy of of the
-  // corresponding function in step-6. The only
-  // difference is the call to
-  // <code>assemble_multigrid</code> that takes
-  // care of forming the matrices on every
-  // level that we need in the multigrid
-  // method.
-  template <int dim>
-  void LaplaceProblem<dim>::run ()
-  {
-    for (unsigned int cycle=0; cycle<15; ++cycle)
-      {
-        pcout << "Cycle " << cycle << ':' << std::endl;
-
-        if (cycle == 0)
-          {
-            GridGenerator::hyper_cube (triangulation);
-
-            triangulation.refine_global (4);
-          }
-        else
-          refine_grid ();
-
-        pcout << "   Number of active cells:       "
-              << triangulation.n_global_active_cells()
-              << std::endl;
-
-        setup_system ();
-
-        pcout << "   Number of degrees of freedom: "
-              << mg_dof_handler.n_dofs()
-              << " (by level: ";
-        for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
-          pcout << mg_dof_handler.n_dofs(level)
-                << (level == triangulation.n_global_levels()-1
-                    ? ")" : ", ");
-        pcout << std::endl;
-
-        assemble_system ();
-        assemble_multigrid ();
-
-        solve ();
-        output_results (cycle);
-      }
-  }
-}
-
-
-// @sect3{The main() function}
-//
-// This is again the same function as
-// in step-6:
-int main (int argc, char *argv[])
-{
-  dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
-
-  try
-    {
-      using namespace dealii;
-      using namespace Step50;
-
-      LaplaceProblem<2> laplace_problem(1/*degree*/);
-      laplace_problem.run ();
-    }
-  catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-                << exc.what() << std::endl
-                << "Aborting!" << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      throw;
-    }
-  catch (...)
-    {
-      std::cerr << std::endl << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-                << "Aborting!" << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      throw;
-    }
-
-  return 0;
-}

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.