]> https://gitweb.dealii.org/ - dealii.git/commitdiff
remove tests/benchmarks/ 2759/head
authorTimo Heister <timo.heister@gmail.com>
Thu, 7 Jul 2016 11:19:28 +0000 (12:19 +0100)
committerTimo Heister <timo.heister@gmail.com>
Thu, 7 Jul 2016 11:19:28 +0000 (12:19 +0100)
they have been located at https://github.com/dealii/performance-
benchmarks for a long time

23 files changed:
tests/benchmarks/baselineplot.sh [deleted file]
tests/benchmarks/bench.sh [deleted file]
tests/benchmarks/benchrev.sh [deleted file]
tests/benchmarks/doplots.sh [deleted file]
tests/benchmarks/gettimes/Makefile [deleted file]
tests/benchmarks/gettimes/get_times.cc [deleted file]
tests/benchmarks/interactive.py [deleted file]
tests/benchmarks/plot.sh [deleted file]
tests/benchmarks/qsub.regression [deleted file]
tests/benchmarks/run.sh [deleted file]
tests/benchmarks/setup.sh [deleted file]
tests/benchmarks/step-22/CMakeLists.txt [deleted file]
tests/benchmarks/step-22/step-22.cc [deleted file]
tests/benchmarks/sync.sh [deleted file]
tests/benchmarks/tablehandler/CMakeLists.txt [deleted file]
tests/benchmarks/tablehandler/table_handler.cc [deleted file]
tests/benchmarks/test_assembly/CMakeLists.txt [deleted file]
tests/benchmarks/test_assembly/step-22.cc [deleted file]
tests/benchmarks/test_hp/CMakeLists.txt [deleted file]
tests/benchmarks/test_hp/step.cc [deleted file]
tests/benchmarks/test_poisson/CMakeLists.txt [deleted file]
tests/benchmarks/test_poisson/poisson.cc [deleted file]
tests/benchmarks/testlist.sh [deleted file]

diff --git a/tests/benchmarks/baselineplot.sh b/tests/benchmarks/baselineplot.sh
deleted file mode 100755 (executable)
index 8ae6bd3..0000000
+++ /dev/null
@@ -1,36 +0,0 @@
-#!/bin/bash
-
-
-source testlist.sh
-
-
-cat <<EOF
-set terminal postscript eps color enh
-set key left top
-set output "baseline.eps"
-#set log y
-set xlabel 'revision'
-set ylabel 'percentage slowdown to baseline'
-set title 'benchmark baselines'
-EOF
-#echo "set terminal x11 persist"
-echo "set xrange [28000:*]"
-echo "set yrange [*:300]"
-echo "plot \\"
-
-for test in $TESTS ; do
-    
-    col=1
-    while read line;
-    do
-       col=`expr $col "+" 1`
-       baseline=`head -n 1 datatable.$test | cut -f $col -d ' '`
-       echo "'datatable.$test' using 1:((\$$col-$baseline)/$baseline*100.0) title '$test - $line' w lp,\\";
-    done < names.$test
-
-done
-
-echo "0.0 w l title 'baseline'"
-
-#echo "pause -1"
-
diff --git a/tests/benchmarks/bench.sh b/tests/benchmarks/bench.sh
deleted file mode 100755 (executable)
index 347f7c0..0000000
+++ /dev/null
@@ -1,28 +0,0 @@
-#!/bin/bash
-
-PREVREVISION="`svn info deal.II | grep Revision | sed s/Revision://`"
-HEADREVISION="`svn info https://svn.dealii.org/trunk/deal.II | grep Revision | sed s/Revision://`"
-
-echo "previous $PREVREVISION"
-echo "HEAD: $HEADREVISION"
-
-while [ $PREVREVISION -lt $HEADREVISION ] ; do
-
-  NEXTREVISION=`expr $PREVREVISION "+" 25`
-  echo "Updating from $PREVREVISION to $NEXTREVISION"
-  cd deal.II
-  svn up -r$NEXTREVISION || exit 1
-  if test -z "`svn diff -r$PREVREVISION:$NEXTREVISION .`" ; then
-      echo "Skipping revision $NEXTREVISION" ;
-      PREVREVISION=$NEXTREVISION
-      cd ..
-      continue ;
-  fi
-
-  cd ..      
-  PREVREVISION=$NEXTREVISION
-  . benchrev.sh 
-
-done
-
-echo "DONE WITH REGRESSION TESTS ON `date`"
diff --git a/tests/benchmarks/benchrev.sh b/tests/benchmarks/benchrev.sh
deleted file mode 100755 (executable)
index 45a963e..0000000
+++ /dev/null
@@ -1,40 +0,0 @@
-#!/bin/bash
-
-source testlist.sh
-
-PREVREVISION="`svn info deal.II | grep Revision | sed s/Revision://`"
-MAKECMD="nice make -j10 install"
-#MAKECMD="nice make -j10 optimized"
-export MAKECMD
-
-echo "testing $PREVREVISION"
-
-  echo "configure"
-  cd build
-  cmake ../deal.II || exit 2
-  echo "compiling" 
-  $MAKECMD || exit 3
-  
-  cd ..
-
-  for test in $TESTS ; do
-      cd $test
-      echo "** working on $test"
-      make clean >/dev/null
-      echo -n "" > temp.txt
-      for a in {1..5}; do
-         echo "*" >> temp.txt
-          make run | grep "|" >> temp.txt
-      done
-      ./../gettimes/gettimes > names.test
-      if [[ -s names.test ]] ; then
-      words=`wc -w names.test | cut -f1 -d' '`
-      if [ "$words" -gt "0" ] ; then
-         cp names.test ../names.$test
-      fi ;
-      rm -rf names.test
-      fi ;
-      ./../gettimes/gettimes $PREVREVISION >>../datatable.$test
-      cd ..      
-  
-  done  
diff --git a/tests/benchmarks/doplots.sh b/tests/benchmarks/doplots.sh
deleted file mode 100755 (executable)
index 7358c04..0000000
+++ /dev/null
@@ -1,20 +0,0 @@
-#!/bin/bash
-
-source testlist.sh
-
-echo "generating images..."
-
-for test in $TESTS ; do
-    LASTREV=`tail -n 1 datatable.$test | cut -f 1 -d ' '`
-    ./plot.sh $test $LASTREV >script
-    gnuplot script
-    rm -rf script
-    convert -density 150 $test.eps $test.png
-done
-
-
-./baselineplot.sh > script
-gnuplot script
-rm -rf script
-convert -density 150 baseline.eps baseline.png
-python interactive.py >index.html
\ No newline at end of file
diff --git a/tests/benchmarks/gettimes/Makefile b/tests/benchmarks/gettimes/Makefile
deleted file mode 100644 (file)
index 3efb000..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-
-gettimes: get_times.cc Makefile
-       @g++ get_times.cc -o gettimes
diff --git a/tests/benchmarks/gettimes/get_times.cc b/tests/benchmarks/gettimes/get_times.cc
deleted file mode 100644 (file)
index 35d4ab8..0000000
+++ /dev/null
@@ -1,166 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2013 - 2014 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.
-//
-// ---------------------------------------------------------------------
-
-
-#include <iostream>
-#include <string>
-#include <sstream>
-#include <fstream>
-#include <vector>
-#include <assert.h>
-#include <stdlib.h>
-
-using namespace std;
-
-int main(int argc, char *argv[])
-{
-
-  ifstream input;
-  stringstream ss;
-  const int IGNORE_LINES = 3; //Number of initial lines to ignore
-  const string DELIM = "*"; //String to divide test data
-
-  input.open("temp.txt");
-
-  vector<string> names;
-  vector<double> times;
-
-  //If no revision number, retrieve column names
-  if (argc <= 1)
-    {
-      string curr_line;
-
-      //Ignore first n lines + delimeter- they don't hold any useful data
-      for (int i = 0; i < IGNORE_LINES + 1; i++)
-        {
-          getline(input, curr_line);
-        }
-
-      while (!input.eof())
-        {
-
-          getline(input, curr_line);
-          if (curr_line == "")
-            continue;
-          if (curr_line == DELIM)
-            break;
-
-          ////cout << "curr line: " << curr_line << endl;
-          //Looking for the string after the '|' and ' '
-          int first_char = 0;
-          if (curr_line[0] == '|')
-            first_char++;
-          if (curr_line[1] == ' ')
-            first_char++;
-
-
-          ////cout << "first char at pos: " << first_char << endl;
-
-          //Find end of string
-          int num_chars = 0;
-          while (curr_line[first_char + num_chars] != '|')
-            num_chars++;
-
-          num_chars--;
-
-          while (curr_line[first_char + num_chars] == ' ')
-            num_chars--;
-
-          names.push_back(curr_line.substr(first_char, num_chars+1));
-        }
-    }
-
-  else //Else, extract execution time from each line
-    {
-
-      int time_index = 0;
-
-      while (!input.eof())
-        {
-
-          string curr_line = "";
-
-          //Read in line
-          getline(input,curr_line);
-
-          //Check for delimeter
-          if (curr_line == DELIM)
-            {
-              ////cout << "Delimeter detected" << endl;
-              time_index = 0;
-              //Ignore first n lines- they don't hold any useful data
-              for (int i = 0; i < IGNORE_LINES; i++)
-                {
-                  string dummy;
-                  getline(input, dummy);
-                  ////cout << "Skipping: " << dummy << endl;
-                }
-              continue;
-            }
-
-          ////cout << "Reading: " << curr_line << endl;
-
-          if (curr_line == "" && DELIM != "")
-            continue;
-
-
-          //Looking for a number that ends with 's'
-          int last_s = curr_line.rfind('s');
-          //In case 's' is not used in the future
-          assert(isdigit(curr_line[last_s - 1])); // Test: s preceded by number
-
-          //Find time string and convert to double
-          int num_start = last_s - 1;
-          while (curr_line[num_start] != ' ')
-            num_start--;
-          string timestr = curr_line.substr(num_start+1, last_s - num_start);
-          double time = (double)atof(timestr.substr(0,timestr.size()-1).c_str());
-
-          assert(times.size() >= time_index);
-          // first addition of times to vector; each loop, times.size() should be one less than time_index
-          if (times.size() == time_index)
-            times.push_back(time);
-          else
-            times[time_index] = (time < times[time_index]) ? time : times[time_index]; //else, determines minimum time and stores it
-          time_index++;
-
-        }
-    }
-
-
-  //Output individual names
-  if (argc <= 1)
-    {
-      for (int i = 0; i < names.size(); i++)
-        {
-          cout << names[i];
-          if (i < names.size() - 1)
-            cout << endl;
-        }
-    }
-
-  //Output individual times
-  if (argc > 1)
-    {
-      cout << argv[1];
-      for (int i = 0; i < times.size(); i++)
-        cout << " " << times[i];
-    }
-
-  cout << endl;
-
-
-  return 0;
-}
diff --git a/tests/benchmarks/interactive.py b/tests/benchmarks/interactive.py
deleted file mode 100755 (executable)
index f687f4b..0000000
+++ /dev/null
@@ -1,155 +0,0 @@
-import textwrap
-import os
-
-begin = \
-"""
-<!DOCTYPE HTML>
-<html>
-       <head>
-               <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
-               <title>deal.II regression timings</title>
-
-               <script type="text/javascript" src="http://ajax.googleapis.com/ajax/libs/jquery/1.8.2/jquery.min.js"></script>
-
-<script src="http://code.highcharts.com/highcharts.js"></script>
-<script src="http://code.highcharts.com/modules/exporting.js"></script>
-
-<script>
-$(function () {
-    var chart;
-    $(document).ready(function() {
-        chart = new Highcharts.Chart({
-            chart: {
-                renderTo: 'container',
-                type: 'line',
-                marginRight: 250,
-                marginBottom: 25,
-                zoomType: 'x'
-            },
-            title: {
-                text: 'regression timings',
-                x: -20 //center
-            },
-            xAxis: {
-            },
-            yAxis: {
-                title: {
-                    text: 'slowdown (%)'
-                },
-                plotLines: [{
-                    value: 0,
-                    width: 1,
-                    color: '#808080'
-                }]
-            },
-            tooltip: {
-                formatter: function() {
-                        return '<b>'+ this.series.name +'</b><br/>'+
-                        this.y +'s' + ', rev ' + this.x;
-                }
-            },
-            legend: {
-                layout: 'vertical',
-                align: 'right',
-                verticalAlign: 'top',
-                x: -10,
-                y: 100,
-                borderWidth: 0
-            },
-
-            series: [
-"""
-
-end = \
-"""
-]
-        });
-setTimeout(function(){
-chart.xAxis[0].setExtremes(28000,chart.xAxis[0].getExtremes().dataMax); 
-},1000);
-    });
-});    
-
-
-</script>
-       </head>
-       <body>
-deal.II performance benchmarks, see 
-<a href="http://www.dealii.org/testsuite.html">http://www.dealii.org/testsuite.html</a><br><br>
-
-<div id="container" style="min-width: 400px; height: 400px; margin: 0 auto"></div>
-
-<img src="baseline.png"/>
-
-<img src="test_assembly.png"/>
-
-<img src="step-22.png"/>
-
-<img src="test_poisson.png"/>
-
-<img src="tablehandler.png"/>
-
-
-       </body>
-</html>
-"""
-
-
-
-list = os.listdir(".")
-
-print begin
-
-first = 1
-for fname in list:
-    if (fname.startswith("names.")):
-       testname = fname[6:]
-       names = open(fname).readlines()
-       data = open("datatable."+testname).readlines()
-       idx = 0
-       for name in names:
-           if first == 1:
-               first = 0
-           else:
-               print ","
-           idx = idx+1
-           print "{ name: '%s - %s', data: [" % (testname,name[:-1])
-           i=0
-#((\$$col-$baseline)/$baseline*100.0)
-           baseline = -1
-           for l in data:
-               if (len(l.strip())<1):
-                   continue
-
-                if (baseline>-1 and baseline<0.5):
-                    continue;
-               
-                def isfloat(x):
-                    try:
-                        float(x)
-                        return True
-                    except ValueError:
-                        return False
-               lnumbers = [float(x) for x in l.split() if isfloat(x)]
-               if len(lnumbers)<=1:
-                   continue;
-
-               if baseline<0 and len(lnumbers)>idx:
-                   baseline=lnumbers[idx]
-
-               if lnumbers[0]<27000:
-                   continue;
-                
-               if (i==1):
-                   print(","),
-               i=1;
-
-                   
-                if len(lnumbers)>idx:
-                    print "[%d,%f]" % (lnumbers[0], (lnumbers[idx]-baseline)/baseline*100.0)
-           print "]}\n"
-       
-
-
-
-print end
diff --git a/tests/benchmarks/plot.sh b/tests/benchmarks/plot.sh
deleted file mode 100755 (executable)
index ec9b452..0000000
+++ /dev/null
@@ -1,28 +0,0 @@
-#!/bin/bash
-
-#launch with the name of test to generate the .eps for
-
-cat <<EOF
-set terminal postscript eps color enh
-set key left bottom
-set output "$1.eps"
-#set log y
-set xlabel 'revision'
-set title 'benchmark $1 - rev $2'
-EOF
-#echo "set terminal x11 persist"
-
-echo "plot \\"
-n=1
-while read line;
-do
-  n=`expr $n "+" 1`
-#  echo "'datatable.$1' using 1:(int(\$$n*10.0+0.5)/10.0) title '$line' w lp,\\";
-  echo "'datatable.$1' using 1:$n title '$line' w lp,\\";
-done < names.$1
-
-# this forces 0.01 to be in the yrange and ends the plot list (trailing comma above)
-echo "0.01 title ''"
-
-#echo "pause -1"
-
diff --git a/tests/benchmarks/qsub.regression b/tests/benchmarks/qsub.regression
deleted file mode 100644 (file)
index 7267a3a..0000000
+++ /dev/null
@@ -1,88 +0,0 @@
-#PBS -l nodes=1:mem64gb
-#PBS -l walltime=4:00:00
-#PBS -q c0541
-#PBS -j oe
-#PBS -N deal.II-regression-tests
-
-MAKECMD="make -j16"
-export MAKECMD
-
-function bdie () {
-        echo "Error: $@"
-        exit 1
-}
-
-if test `hostname` != "c0541" ; then
-  echo "Wrong machine!"
-  echo "This is `hostname`..."
-  exit 1
-fi
-
-
-
-echo "STARTING REGRESSION TESTS ON `date`"
-
-cd /node/bangerth/regression-test-do/deal.II/ || bdie "$LINENO"
-
-PREVREVISION="`svn info . | grep Revision | sed s/Revision://`"
-HEADREVISION="`svn info http://www.dealii.org/svn/dealii | grep Revision | sed s/Revision://`"
-
-if test "$PREVREVISION" = "$HEADREVISION" ; then
-  echo "$PREVREVISION already handled"
-else
-  NEXTREVISION=`expr $PREVREVISION "+" 1`
-  echo "Updating from $PREVREVISION $NEXTREVISION"
-
-  for dir in /node/bangerth/regression-test-do/deal.II \
-             /node/bangerth/regression-test-do/branch_higher_derivatives/deal.II \
-            /node/bangerth/regression-test-do/branch_merge_mg_into_dof_handler/deal.II \
-            /node/bangerth/regression-test-do/branch_component_mask/deal.II ; do
-    cd $dir ;
-    svn up -r$NEXTREVISION   || bdie "$LINENO"
-    svn up tests -r$NEXTREVISION || bdie "$LINENO"
-
-    # see if anything changed between the previous revision
-    # and the current one. if not, then simply skip this branch
-    if test -z "`svn diff -r$PREVREVISION:$NEXTREVISION . tests`" ; then
-      echo "Skipping revision $NEXTREVISION on branch $dir" ;
-      continue ;
-    else
-      echo "Doing revision $NEXTREVISION on branch $dir because of these diffs:" ;
-      echo svn diff -r$PREVREVISION:$NEXTREVISION . tests ;
-      svn diff -r$PREVREVISION:$NEXTREVISION . tests ;
-    fi
-
-    rm -f source/Makefile.dep
-    # $MAKECMD clean distclean || bdie "$LINENO"
-
-    # use libstdc++ debug mode for our tests
-    #  export CXXFLAGS=-D_GLIBCXX_DEBUG
-    ./reconfigure            || bdie "$LINENO"
-    $MAKECMD debug           || bdie "$LINENO"
-    cd tests                 || bdie "$LINENO"
-    $MAKECMD clean distclean
-    $MAKECMD report+mail     || bdie "$LINENO"
-    touch sent-mail-for-revision-$NEXTREVISION
-
-    if test $dir = /node/bangerth/regression-test-do/deal.II ; then
-      cd ../../projects/aspect
-      svn up
-      $MAKECMD clean
-      $MAKECMD
-      cd tests
-      $MAKECMD report+mail     || bdie "$LINENO"
-    fi
-  done
-
-
-  cd
-  if test "$NEXTREVISION" = "$HEADREVISION" ; then
-    : ;
-  else
-    echo "Chaining next revision after $NEXTREVISION (current is $HEADREVISION)"
-    qsub qsub.regression
-  fi
-fi
-
-echo "DONE WITH REGRESSION TESTS ON `date`"
-
diff --git a/tests/benchmarks/run.sh b/tests/benchmarks/run.sh
deleted file mode 100755 (executable)
index 939c670..0000000
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/bin/bash
-
-
-while :
-do
-   ./bench.sh
-
-   ./doplots.sh
-   
-   ./sync.sh
-
-   sleep 3600
-done
\ No newline at end of file
diff --git a/tests/benchmarks/setup.sh b/tests/benchmarks/setup.sh
deleted file mode 100755 (executable)
index a8224ef..0000000
+++ /dev/null
@@ -1,28 +0,0 @@
-#!/bin/bash
-
-# set the version we want to start with. the first revision that
-# can be used with this script is r24500
-REV=30125
-
-rm -rf deal.II
-svn co -r $REV https://svn.dealii.org/trunk/deal.II
-
-mkdir build
-cd build
-cmake -DCMAKE_BUILD_TYPE=Release -DDEAL_II_WITH_THREADS=OFF -DCMAKE_INSTALL_PREFIX=`pwd`/../installed ../deal.II
-make install -j 10
-cd ..
-
-cd gettimes
-make
-cd ..
-
-source testlist.sh
-for test in $TESTS ; do
-      cd $test
-      echo "** cmake for $test:"
-      cmake -DDEAL_II_DIR=`pwd`/../installed
-      make release
-      cd ..
-done
-
diff --git a/tests/benchmarks/step-22/CMakeLists.txt b/tests/benchmarks/step-22/CMakeLists.txt
deleted file mode 100644 (file)
index f261836..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-##
-#  CMake script for the step-1 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "step-22")
-
-# Declare all source files the target consists of:
-SET(TARGET_SRC
-  ${TARGET}.cc
-  # You can specify additional files here!
-  )
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.0 QUIET
-  HINTS
-    ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR}
-  #
-  # If the deal.II library cannot be found (because it is not installed at a
-  # default location or your project resides at an uncommon place), you
-  # can specify additional hints for search paths here, e.g.
-  # "$ENV{HOME}/workspace/deal.II"
-  )
-
-IF (NOT ${deal.II_FOUND})
-   MESSAGE(FATAL_ERROR
-           "\n\n"
-          " *** Could not locate deal.II. *** "
-          "\n\n"
-           " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag 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/tests/benchmarks/step-22/step-22.cc b/tests/benchmarks/step-22/step-22.cc
deleted file mode 100644 (file)
index 7753597..0000000
+++ /dev/null
@@ -1,1386 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2008 - 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.
-//
-// ---------------------------------------------------------------------
-
-/*
- * step-22.cc
- */
-
-// @sect3{Include files}
-
-// As usual, we start by including
-// some well-known files:
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/constraint_matrix.h>
-
-#include <deal.II/grid/tria.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/tria_boundary_lib.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/grid_refinement.h>
-
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_renumbering.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_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 <deal.II/numerics/error_estimator.h>
-
-// Then we need to include the header file
-// for the sparse direct solver UMFPACK:
-#include <deal.II/lac/sparse_direct.h>
-#include <deal.II/base/timer.h>
-
-// This includes the library for the
-// incomplete LU factorization that will
-// be used as a preconditioner in 3D:
-#include <deal.II/lac/sparse_ilu.h>
-
-// This is C++:
-#include <fstream>
-#include <sstream>
-
-// As in all programs, the namespace dealii
-// is included:
-namespace Step22
-{
-  using namespace dealii;
-
-  // @sect3{Defining the inner preconditioner type}
-
-  // As explained in the introduction, we are
-  // going to use different preconditioners for
-  // two and three space dimensions,
-  // respectively. We distinguish between
-  // them by the use of the spatial dimension
-  // as a template parameter. See step-4 for
-  // details on templates. We are not going to
-  // create any preconditioner object here, all
-  // we do is to create class that holds a
-  // local typedef determining the
-  // preconditioner class so we can write our
-  // program in a dimension-independent way.
-  template <int dim>
-  struct InnerPreconditioner;
-
-  // In 2D, we are going to use a sparse direct
-  // solver as preconditioner:
-  template <>
-  struct InnerPreconditioner<2>
-  {
-    typedef SparseILU<double> type;
-//      typedef SparseDirectUMFPACK type;
-  };
-
-  // And the ILU preconditioning in 3D, called
-  // by SparseILU:
-  template <>
-  struct InnerPreconditioner<3>
-  {
-    typedef SparseILU<double> type;
-  };
-
-
-  // @sect3{The <code>StokesProblem</code> class template}
-
-  // This is an adaptation of step-20, so the
-  // main class and the data types are the
-  // same as used there. In this example we
-  // also use adaptive grid refinement, which
-  // is handled in analogy to
-  // step-6. According to the discussion in
-  // the introduction, we are also going to
-  // use the ConstraintMatrix for
-  // implementing Dirichlet boundary
-  // conditions. Hence, we change the name
-  // <code>hanging_node_constraints</code>
-  // into <code>constraints</code>.
-  template <int dim>
-  class StokesProblem
-  {
-  public:
-    StokesProblem (const unsigned int degree);
-    void run ();
-
-  private:
-    void setup_dofs ();
-    void assemble_system ();
-    void solve ();
-    void output_results (const unsigned int refinement_cycle) const;
-    void refine_mesh ();
-
-    const unsigned int   degree;
-
-    Triangulation<dim>   triangulation;
-    FESystem<dim>        fe;
-    DoFHandler<dim>      dof_handler;
-
-    ConstraintMatrix     constraints;
-
-    BlockSparsityPattern      sparsity_pattern;
-    BlockSparseMatrix<double> system_matrix;
-
-    BlockVector<double> solution;
-    BlockVector<double> system_rhs;
-
-    // This one is new: We shall use a
-    // so-called shared pointer structure to
-    // access the preconditioner. Shared
-    // pointers are essentially just a
-    // convenient form of pointers. Several
-    // shared pointers can point to the same
-    // object (just like regular pointers),
-    // but when the last shared pointer
-    // object to point to a preconditioner
-    // object is deleted (for example if a
-    // shared pointer object goes out of
-    // scope, if the class of which it is a
-    // member is destroyed, or if the pointer
-    // is assigned a different preconditioner
-    // object) then the preconditioner object
-    // pointed to is also destroyed. This
-    // ensures that we don't have to manually
-    // track in how many places a
-    // preconditioner object is still
-    // referenced, it can never create a
-    // memory leak, and can never produce a
-    // dangling pointer to an already
-    // destroyed object:
-    std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
-
-    TimerOutput timer;
-  };
-
-  // @sect3{Boundary values and right hand side}
-
-  // As in step-20 and most other
-  // example programs, the next task is
-  // to define the data for the PDE:
-  // For the Stokes problem, we are
-  // going to use natural boundary
-  // values on parts of the boundary
-  // (i.e. homogenous Neumann-type) for
-  // which we won't have to do anything
-  // special (the homogeneity implies
-  // that the corresponding terms in
-  // the weak form are simply zero),
-  // and boundary conditions on the
-  // velocity (Dirichlet-type) on the
-  // rest of the boundary, as described
-  // in the introduction.
-  //
-  // In order to enforce the Dirichlet
-  // boundary values on the velocity,
-  // we will use the
-  // VectorTools::interpolate_boundary_values
-  // function as usual which requires
-  // us to write a function object with
-  // as many components as the finite
-  // element has. In other words, we
-  // have to define the function on the
-  // $(u,p)$-space, but we are going to
-  // filter out the pressure component
-  // when interpolating the boundary
-  // values.
-
-  // The following function object is a
-  // representation of the boundary
-  // values described in the
-  // introduction:
-  template <int dim>
-  class BoundaryValues : public Function<dim>
-  {
-  public:
-    BoundaryValues () : Function<dim>(dim+1) {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
-  };
-
-
-  template <int dim>
-  double
-  BoundaryValues<dim>::value (const Point<dim>  &p,
-                              const unsigned int component) const
-  {
-    Assert (component < this->n_components,
-            ExcIndexRange (component, 0, this->n_components));
-
-    if (component == 0)
-      return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
-    return 0;
-  }
-
-
-  template <int dim>
-  void
-  BoundaryValues<dim>::vector_value (const Point<dim> &p,
-                                     Vector<double>   &values) const
-  {
-    for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = BoundaryValues<dim>::value (p, c);
-  }
-
-
-
-  // We implement similar functions for
-  // the right hand side which for the
-  // current example is simply zero:
-  template <int dim>
-  class RightHandSide : public Function<dim>
-  {
-  public:
-    RightHandSide () : Function<dim>(dim+1) {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
-
-  };
-
-
-  template <int dim>
-  double
-  RightHandSide<dim>::value (const Point<dim>  &/*p*/,
-                             const unsigned int /*component*/) const
-  {
-    return 0;
-  }
-
-
-  template <int dim>
-  void
-  RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                    Vector<double>   &values) const
-  {
-    for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = RightHandSide<dim>::value (p, c);
-  }
-
-
-  // @sect3{Linear solvers and preconditioners}
-
-  // The linear solvers and preconditioners are
-  // discussed extensively in the
-  // introduction. Here, we create the
-  // respective objects that will be used.
-
-  // @sect4{The <code>InverseMatrix</code> class template}
-
-  // The <code>InverseMatrix</code>
-  // class represents the data
-  // structure for an inverse
-  // matrix. It is derived from the one
-  // in step-20. The only difference is
-  // that we now do include a
-  // preconditioner to the matrix since
-  // we will apply this class to
-  // different kinds of matrices that
-  // will require different
-  // preconditioners (in step-20 we did
-  // not use a preconditioner in this
-  // class at all). The types of matrix
-  // and preconditioner are passed to
-  // this class via template
-  // parameters, and matrix and
-  // preconditioner objects of these
-  // types will then be passed to the
-  // constructor when an
-  // <code>InverseMatrix</code> object
-  // is created. The member function
-  // <code>vmult</code> is, as in
-  // step-20, a multiplication with a
-  // vector, obtained by solving a
-  // linear system:
-  template <class Matrix, class Preconditioner>
-  class InverseMatrix : public Subscriptor
-  {
-  public:
-    InverseMatrix (const Matrix         &m,
-                   const Preconditioner &preconditioner);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const Matrix> matrix;
-    const SmartPointer<const Preconditioner> preconditioner;
-  };
-
-
-  template <class Matrix, class Preconditioner>
-  InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
-                                                       const Preconditioner &preconditioner)
-    :
-    matrix (&m),
-    preconditioner (&preconditioner)
-  {}
-
-
-  // This is the implementation of the
-  // <code>vmult</code> function.
-
-  // In this class we use a rather large
-  // tolerance for the solver control. The
-  // reason for this is that the function is
-  // used very frequently, and hence, any
-  // additional effort to make the residual
-  // in the CG solve smaller makes the
-  // solution more expensive. Note that we do
-  // not only use this class as a
-  // preconditioner for the Schur complement,
-  // but also when forming the inverse of the
-  // Laplace matrix &ndash; which is hence
-  // directly responsible for the accuracy of
-  // the solution itself, so we can't choose
-  // a too large tolerance, either.
-  template <class Matrix, class Preconditioner>
-  void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
-                                                    const Vector<double> &src) const
-  {
-    SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
-    SolverCG<>    cg (solver_control);
-
-    dst = 0;
-
-    cg.solve (*matrix, dst, src, *preconditioner);
-  }
-
-
-  // @sect4{The <code>SchurComplement</code> class template}
-
-  // This class implements the Schur complement
-  // discussed in the introduction.  It is in
-  // analogy to step-20.  Though, we now call
-  // it with a template parameter
-  // <code>Preconditioner</code> in order to
-  // access that when specifying the respective
-  // type of the inverse matrix class. As a
-  // consequence of the definition above, the
-  // declaration <code>InverseMatrix</code> now
-  // contains the second template parameter
-  // for a preconditioner class as above, which
-  // affects the <code>SmartPointer</code>
-  // object <code>m_inverse</code> as well.
-  template <class Preconditioner>
-  class SchurComplement : public Subscriptor
-  {
-  public:
-    SchurComplement (const BlockSparseMatrix<double> &system_matrix,
-                     const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
-
-    mutable Vector<double> tmp1, tmp2;
-  };
-
-
-
-  template <class Preconditioner>
-  SchurComplement<Preconditioner>::
-  SchurComplement (const BlockSparseMatrix<double> &system_matrix,
-                   const InverseMatrix<SparseMatrix<double>,Preconditioner> &A_inverse)
-    :
-    system_matrix (&system_matrix),
-    A_inverse (&A_inverse),
-    tmp1 (system_matrix.block(0,0).m()),
-    tmp2 (system_matrix.block(0,0).m())
-  {}
-
-
-  template <class Preconditioner>
-  void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
-                                               const Vector<double> &src) const
-  {
-    system_matrix->block(0,1).vmult (tmp1, src);
-    A_inverse->vmult (tmp2, tmp1);
-    system_matrix->block(1,0).vmult (dst, tmp2);
-  }
-
-
-  // @sect3{StokesProblem class implementation}
-
-  // @sect4{StokesProblem::StokesProblem}
-
-  // The constructor of this class
-  // looks very similar to the one of
-  // step-20. The constructor
-  // initializes the variables for the
-  // polynomial degree, triangulation,
-  // finite element system and the dof
-  // handler. The underlying polynomial
-  // functions are of order
-  // <code>degree+1</code> for the
-  // vector-valued velocity components
-  // and of order <code>degree</code>
-  // for the pressure.  This gives the
-  // LBB-stable element pair
-  // $Q_{degree+1}^d\times Q_{degree}$,
-  // often referred to as the
-  // Taylor-Hood element.
-  //
-  // Note that we initialize the triangulation
-  // with a MeshSmoothing argument, which
-  // ensures that the refinement of cells is
-  // done in a way that the approximation of
-  // the PDE solution remains well-behaved
-  // (problems arise if grids are too
-  // unstructered), see the documentation of
-  // <code>Triangulation::MeshSmoothing</code>
-  // for details.
-  template <int dim>
-  StokesProblem<dim>::StokesProblem (const unsigned int degree)
-    :
-    degree (degree),
-    triangulation (Triangulation<dim>::maximum_smoothing),
-    fe (FE_Q<dim>(degree+1), dim,
-        FE_Q<dim>(degree), 1),
-    dof_handler (triangulation),
-    timer (std::cout, TimerOutput::summary, TimerOutput::cpu_times)
-  {}
-
-
-  // @sect4{StokesProblem::setup_dofs}
-
-  // Given a mesh, this function
-  // associates the degrees of freedom
-  // with it and creates the
-  // corresponding matrices and
-  // vectors. At the beginning it also
-  // releases the pointer to the
-  // preconditioner object (if the
-  // shared pointer pointed at anything
-  // at all at this point) since it
-  // will definitely not be needed any
-  // more after this point and will
-  // have to be re-computed after
-  // assembling the matrix, and unties
-  // the sparse matrix from its
-  // sparsity pattern object.
-  //
-  // We then proceed with distributing
-  // degrees of freedom and renumbering
-  // them: In order to make the ILU
-  // preconditioner (in 3D) work
-  // efficiently, it is important to
-  // enumerate the degrees of freedom
-  // in such a way that it reduces the
-  // bandwidth of the matrix, or maybe
-  // more importantly: in such a way
-  // that the ILU is as close as
-  // possible to a real LU
-  // decomposition. On the other hand,
-  // we need to preserve the block
-  // structure of velocity and pressure
-  // already seen in in step-20 and
-  // step-21. This is done in two
-  // steps: First, all dofs are
-  // renumbered to improve the ILU and
-  // then we renumber once again by
-  // components. Since
-  // <code>DoFRenumbering::component_wise</code>
-  // does not touch the renumbering
-  // within the individual blocks, the
-  // basic renumbering from the first
-  // step remains. As for how the
-  // renumber degrees of freedom to
-  // improve the ILU: deal.II has a
-  // number of algorithms that attempt
-  // to find orderings to improve ILUs,
-  // or reduce the bandwidth of
-  // matrices, or optimize some other
-  // aspect. The DoFRenumbering
-  // namespace shows a comparison of
-  // the results we obtain with several
-  // of these algorithms based on the
-  // testcase discussed here in this
-  // tutorial program. Here, we will
-  // use the traditional Cuthill-McKee
-  // algorithm already used in some of
-  // the previous tutorial programs.
-  // In the
-  // <a href="#improved-ilu">section on improved ILU</a>
-  // we're going to discuss this issue
-  // in more detail.
-
-  // There is one more change compared
-  // to previous tutorial programs:
-  // There is no reason in sorting the
-  // <code>dim</code> velocity
-  // components individually. In fact,
-  // rather than first enumerating all
-  // $x$-velocities, then all
-  // $y$-velocities, etc, we would like
-  // to keep all velocities at the same
-  // location together and only
-  // separate between velocities (all
-  // components) and pressures. By
-  // default, this is not what the
-  // DoFRenumbering::component_wise
-  // function does: it treats each
-  // vector component separately; what
-  // we have to do is group several
-  // components into "blocks" and pass
-  // this block structure to that
-  // function. Consequently, we
-  // allocate a vector
-  // <code>block_component</code> with
-  // as many elements as there are
-  // components and describe all
-  // velocity components to correspond
-  // to block 0, while the pressure
-  // component will form block 1:
-  template <int dim>
-  void StokesProblem<dim>::setup_dofs ()
-  {
-    A_preconditioner.reset ();
-    system_matrix.clear ();
-
-    dof_handler.distribute_dofs (fe);
-    DoFRenumbering::Cuthill_McKee (dof_handler);
-
-    std::vector<unsigned int> block_component (dim+1,0);
-    block_component[dim] = 1;
-    DoFRenumbering::component_wise (dof_handler, block_component);
-
-    // Now comes the implementation of
-    // Dirichlet boundary conditions, which
-    // should be evident after the discussion
-    // in the introduction. All that changed is
-    // that the function already appears in the
-    // setup functions, whereas we were used to
-    // see it in some assembly routine. Further
-    // down below where we set up the mesh, we
-    // will associate the top boundary where we
-    // impose Dirichlet boundary conditions
-    // with boundary indicator 1.  We will have
-    // to pass this boundary indicator as
-    // second argument to the function below
-    // interpolating boundary values.  There is
-    // one more thing, though.  The function
-    // describing the Dirichlet conditions was
-    // defined for all components, both
-    // velocity and pressure. However, the
-    // Dirichlet conditions are to be set for
-    // the velocity only.  To this end, we use
-    // a <code>component_mask</code> that
-    // filters out the pressure component, so
-    // that the condensation is performed on
-    // velocity degrees of freedom only. Since
-    // we use adaptively refined grids the
-    // constraint matrix needs to be first
-    // filled with hanging node constraints
-    // generated from the DoF handler. Note the
-    // order of the two functions &mdash; we
-    // first compute the hanging node
-    // constraints, and then insert the
-    // boundary values into the constraint
-    // matrix. This makes sure that we respect
-    // H<sup>1</sup> conformity on boundaries
-    // with hanging nodes (in three space
-    // dimensions), where the hanging node
-    // needs to dominate the Dirichlet boundary
-    // values.
-    {
-      constraints.clear ();
-      std::vector<bool> component_mask (dim+1, true);
-      component_mask[dim] = false;
-      DoFTools::make_hanging_node_constraints (dof_handler,
-                                               constraints);
-      VectorTools::interpolate_boundary_values (dof_handler,
-                                                1,
-                                                BoundaryValues<dim>(),
-                                                constraints,
-                                                component_mask);
-    }
-
-    constraints.close ();
-
-    // In analogy to step-20, we count the dofs
-    // in the individual components.  We could
-    // do this in the same way as there, but we
-    // want to operate on the block structure
-    // we used already for the renumbering: The
-    // function
-    // <code>DoFTools::count_dofs_per_block</code>
-    // does the same as
-    // <code>DoFTools::count_dofs_per_component</code>,
-    // but now grouped as velocity and pressure
-    // block via <code>block_component</code>.
-    std::vector<unsigned int> dofs_per_block (2);
-    DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
-    const unsigned int n_u = dofs_per_block[0],
-                       n_p = dofs_per_block[1];
-
-    std::cout << "   Number of active cells: "
-              << triangulation.n_active_cells()
-              << std::endl
-              << "   Number of degrees of freedom: "
-              << dof_handler.n_dofs()
-              << " (" << n_u << '+' << n_p << ')'
-              << std::endl;
-
-    // The next task is to allocate a
-    // sparsity pattern for the system matrix
-    // we will create. We could do this in
-    // the same way as in step-20,
-    // i.e. directly build an object of type
-    // SparsityPattern through
-    // DoFTools::make_sparsity_pattern. However,
-    // there is a major reason not to do so:
-    // In 3D, the function
-    // DoFTools::max_couplings_between_dofs
-    // yields a conservative but rather large
-    // number for the coupling between the
-    // individual dofs, so that the memory
-    // initially provided for the creation of
-    // the sparsity pattern of the matrix is
-    // far too much -- so much actually that
-    // the initial sparsity pattern won't
-    // even fit into the physical memory of
-    // most systems already for
-    // moderately-sized 3D problems, see also
-    // the discussion in step-18.  Instead,
-    // we first build a temporary object that
-    // uses a different data structure that
-    // doesn't require allocating more memory
-    // than necessary but isn't suitable for
-    // use as a basis of SparseMatrix or
-    // BlockSparseMatrix objects; in a second
-    // step we then copy this object into an
-    // object of BlockSparsityPattern. This
-    // is entirely analgous to what we
-    // already did in step-11 and step-18.
-    //
-    // There is one snag again here, though:
-    // it turns out that using the
-    // CompressedSparsityPattern (or the
-    // block version
-    // BlockCompressedSparsityPattern we
-    // would use here) has a bottleneck that
-    // makes the algorithm to build the
-    // sparsity pattern be quadratic in the
-    // number of degrees of freedom. This
-    // doesn't become noticeable until we get
-    // well into the range of several 100,000
-    // degrees of freedom, but eventually
-    // dominates the setup of the linear
-    // system when we get to more than a
-    // million degrees of freedom. This is
-    // due to the data structures used in the
-    // CompressedSparsityPattern class,
-    // nothing that can easily be
-    // changed. Fortunately, there is an easy
-    // solution: the
-    // CompressedSimpleSparsityPattern class
-    // (and its block variant
-    // BlockCompressedSimpleSparsityPattern)
-    // has exactly the same interface, uses a
-    // different %internal data structure and
-    // is linear in the number of degrees of
-    // freedom and therefore much more
-    // efficient for large problems. As
-    // another alternative, we could also
-    // have chosen the class
-    // BlockCompressedSetSparsityPattern that
-    // uses yet another strategy for %internal
-    // memory management. Though, that class
-    // turns out to be more memory-demanding
-    // than
-    // BlockCompressedSimpleSparsityPattern
-    // for this example.
-    //
-    // Consequently, this is the class that
-    // we will use for our intermediate
-    // sparsity representation. All this is
-    // done inside a new scope, which means
-    // that the memory of <code>csp</code>
-    // will be released once the information
-    // has been copied to
-    // <code>sparsity_pattern</code>.
-    {
-      BlockCompressedSimpleSparsityPattern csp (2,2);
-
-      csp.block(0,0).reinit (n_u, n_u);
-      csp.block(1,0).reinit (n_p, n_u);
-      csp.block(0,1).reinit (n_u, n_p);
-      csp.block(1,1).reinit (n_p, n_p);
-
-      csp.collect_sizes();
-
-      DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
-      sparsity_pattern.copy_from (csp);
-    }
-
-    // Finally, the system matrix,
-    // solution and right hand side are
-    // created from the block
-    // structure as in step-20:
-    system_matrix.reinit (sparsity_pattern);
-
-    solution.reinit (2);
-    solution.block(0).reinit (n_u);
-    solution.block(1).reinit (n_p);
-    solution.collect_sizes ();
-
-    system_rhs.reinit (2);
-    system_rhs.block(0).reinit (n_u);
-    system_rhs.block(1).reinit (n_p);
-    system_rhs.collect_sizes ();
-  }
-
-
-  // @sect4{StokesProblem::assemble_system}
-
-  // The assembly process follows the
-  // discussion in step-20 and in the
-  // introduction. We use the well-known
-  // abbreviations for the data structures
-  // that hold the local matrix, right
-  // hand side, and global
-  // numbering of the degrees of freedom
-  // for the present cell.
-  template <int dim>
-  void StokesProblem<dim>::assemble_system ()
-  {
-    system_matrix=0;
-    system_rhs=0;
-
-    QGauss<dim>   quadrature_formula(degree+2);
-
-    FEValues<dim> fe_values (fe, quadrature_formula,
-                             update_values    |
-                             update_quadrature_points  |
-                             update_JxW_values |
-                             update_gradients);
-
-    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
-
-    const unsigned int   n_q_points      = quadrature_formula.size();
-
-    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>       local_rhs (dofs_per_cell);
-
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    const RightHandSide<dim>          right_hand_side;
-    std::vector<Vector<double> >      rhs_values (n_q_points,
-                                                  Vector<double>(dim+1));
-
-    // Next, we need two objects that work as
-    // extractors for the FEValues
-    // object. Their use is explained in detail
-    // in the report on @ref vector_valued :
-    const FEValuesExtractors::Vector velocities (0);
-    const FEValuesExtractors::Scalar pressure (dim);
-
-    // As an extension over step-20 and
-    // step-21, we include a few
-    // optimizations that make assembly
-    // much faster for this particular
-    // problem.  The improvements are
-    // based on the observation that we
-    // do a few calculations too many
-    // times when we do as in step-20:
-    // The symmetric gradient actually
-    // has <code>dofs_per_cell</code>
-    // different values per quadrature
-    // point, but we extract it
-    // <code>dofs_per_cell*dofs_per_cell</code>
-    // times from the FEValues object -
-    // for both the loop over
-    // <code>i</code> and the inner
-    // loop over <code>j</code>. In 3d,
-    // that means evaluating it
-    // $89^2=7921$ instead of $89$
-    // times, a not insignificant
-    // difference.
-    //
-    // So what we're
-    // going to do here is to avoid
-    // such repeated calculations by
-    // getting a vector of rank-2
-    // tensors (and similarly for
-    // the divergence and the basis
-    // function value on pressure)
-    // at the quadrature point prior
-    // to starting the loop over the
-    // dofs on the cell. First, we
-    // create the respective objects
-    // that will hold these
-    // values. Then, we start the
-    // loop over all cells and the loop
-    // over the quadrature points,
-    // where we first extract these
-    // values. There is one more
-    // optimization we implement here:
-    // the local matrix (as well as
-    // the global one) is going to
-    // be symmetric, since all
-    // the operations involved are
-    // symmetric with respect to $i$
-    // and $j$. This is implemented by
-    // simply running the inner loop
-    // not to <code>dofs_per_cell</code>,
-    // but only up to <code>i</code>,
-    // the index of the outer loop.
-    std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
-    std::vector<double>                  div_phi_u   (dofs_per_cell);
-    std::vector<double>                  phi_p       (dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-    for (; cell!=endc; ++cell)
-      {
-        fe_values.reinit (cell);
-        local_matrix = 0;
-        local_rhs = 0;
-
-        right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
-                                          rhs_values);
-
-        for (unsigned int q=0; q<n_q_points; ++q)
-          {
-            for (unsigned int k=0; k<dofs_per_cell; ++k)
-              {
-                symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
-                div_phi_u[k]     = fe_values[velocities].divergence (k, q);
-                phi_p[k]         = fe_values[pressure].value (k, q);
-              }
-
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              {
-                for (unsigned int j=0; j<=i; ++j)
-                  {
-                    local_matrix(i,j) += (symgrad_phi_u[i] * symgrad_phi_u[j]
-                                          - div_phi_u[i] * phi_p[j]
-                                          - phi_p[i] * div_phi_u[j]
-                                          + phi_p[i] * phi_p[j])
-                                         * fe_values.JxW(q);
-
-                  }
-
-                const unsigned int component_i =
-                  fe.system_to_component_index(i).first;
-                local_rhs(i) += fe_values.shape_value(i,q) *
-                                rhs_values[q](component_i) *
-                                fe_values.JxW(q);
-              }
-          }
-
-        // Note that in the above computation
-        // of the local matrix contribution
-        // we added the term <code> phi_p[i] *
-        // phi_p[j] </code>, yielding a
-        // pressure mass matrix in the
-        // $(1,1)$ block of the matrix as
-        // discussed in the
-        // introduction. That this term only
-        // ends up in the $(1,1)$ block stems
-        // from the fact that both of the
-        // factors in <code>phi_p[i] *
-        // phi_p[j]</code> are only non-zero
-        // when all the other terms vanish
-        // (and the other way around).
-        //
-        // Note also that operator* is
-        // overloaded for symmetric
-        // tensors, yielding the scalar
-        // product between the two
-        // tensors in the first line of
-        // the local matrix
-        // contribution.
-
-        // Before we can write the local data
-        // into the global matrix (and
-        // simultaneously use the
-        // ConstraintMatrix object to apply
-        // Dirichlet boundary conditions and
-        // eliminate hanging node
-        // constraints, as we discussed in
-        // the introduction), we have to be
-        // careful about one thing,
-        // though. We have only build up half
-        // of the local matrix because of
-        // symmetry, but we're going to save
-        // the full system matrix in order to
-        // use the standard functions for
-        // solution. This is done by flipping
-        // the indices in case we are
-        // pointing into the empty part of
-        // the local matrix.
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          for (unsigned int j=i+1; j<dofs_per_cell; ++j)
-            local_matrix(i,j) = local_matrix(j,i);
-
-        cell->get_dof_indices (local_dof_indices);
-        constraints.distribute_local_to_global (local_matrix, local_rhs,
-                                                local_dof_indices,
-                                                system_matrix, system_rhs);
-      }
-
-    // Before we're going to solve this
-    // linear system, we generate a
-    // preconditioner for the
-    // velocity-velocity matrix, i.e.,
-    // <code>block(0,0)</code> in the
-    // system matrix. As mentioned
-    // above, this depends on the
-    // spatial dimension. Since the two
-    // classes described by the
-    // <code>InnerPreconditioner::type</code>
-    // typedef have the same interface,
-    // we do not have to do anything
-    // different whether we want to use
-    // a sparse direct solver or an
-    // ILU:
-    std::cout << "   Computing preconditioner..." << std::endl << std::flush;
-
-    A_preconditioner
-      = std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
-    A_preconditioner->initialize (system_matrix.block(0,0),
-                                  typename InnerPreconditioner<dim>::type::AdditionalData());
-
-  }
-
-
-
-  // @sect4{StokesProblem::solve}
-
-  // After the discussion in the introduction
-  // and the definition of the respective
-  // classes above, the implementation of the
-  // <code>solve</code> function is rather
-  // straigt-forward and done in a similar way
-  // as in step-20. To start with, we need an
-  // object of the <code>InverseMatrix</code>
-  // class that represents the inverse of the
-  // matrix A. As described in the
-  // introduction, the inverse is generated
-  // with the help of an inner preconditioner
-  // of type
-  // <code>InnerPreconditioner::type</code>.
-  template <int dim>
-  void StokesProblem<dim>::solve ()
-  {
-    const InverseMatrix<SparseMatrix<double>,
-          typename InnerPreconditioner<dim>::type>
-          A_inverse (system_matrix.block(0,0), *A_preconditioner);
-    Vector<double> tmp (solution.block(0).size());
-
-    // This is as in step-20. We generate the
-    // right hand side $B A^{-1} F - G$ for the
-    // Schur complement and an object that
-    // represents the respective linear
-    // operation $B A^{-1} B^T$, now with a
-    // template parameter indicating the
-    // preconditioner - in accordance with the
-    // definition of the class.
-    {
-      Vector<double> schur_rhs (solution.block(1).size());
-      A_inverse.vmult (tmp, system_rhs.block(0));
-      system_matrix.block(1,0).vmult (schur_rhs, tmp);
-      schur_rhs -= system_rhs.block(1);
-
-      SchurComplement<typename InnerPreconditioner<dim>::type>
-      schur_complement (system_matrix, A_inverse);
-
-      // The usual control structures for
-      // the solver call are created...
-      SolverControl solver_control (solution.block(1).size(),
-                                    1e-6*schur_rhs.l2_norm());
-      SolverCG<>    cg (solver_control);
-
-      // Now to the preconditioner to the
-      // Schur complement. As explained in
-      // the introduction, the
-      // preconditioning is done by a mass
-      // matrix in the pressure variable.  It
-      // is stored in the $(1,1)$ block of
-      // the system matrix (that is not used
-      // anywhere else but in
-      // preconditioning).
-      //
-      // Actually, the solver needs to have
-      // the preconditioner in the form
-      // $P^{-1}$, so we need to create an
-      // inverse operation. Once again, we
-      // use an object of the class
-      // <code>InverseMatrix</code>, which
-      // implements the <code>vmult</code>
-      // operation that is needed by the
-      // solver.  In this case, we have to
-      // invert the pressure mass matrix. As
-      // it already turned out in earlier
-      // tutorial programs, the inversion of
-      // a mass matrix is a rather cheap and
-      // straight-forward operation (compared
-      // to, e.g., a Laplace matrix). The CG
-      // method with ILU preconditioning
-      // converges in 5-10 steps,
-      // independently on the mesh size.
-      // This is precisely what we do here:
-      // We choose another ILU preconditioner
-      // and take it along to the
-      // InverseMatrix object via the
-      // corresponding template parameter.  A
-      // CG solver is then called within the
-      // vmult operation of the inverse
-      // matrix.
-      //
-      // An alternative that is cheaper to
-      // build, but needs more iterations
-      // afterwards, would be to choose a
-      // SSOR preconditioner with factor
-      // 1.2. It needs about twice the number
-      // of iterations, but the costs for its
-      // generation are almost neglible.
-      SparseILU<double> preconditioner;
-      preconditioner.initialize (system_matrix.block(1,1),
-                                 SparseILU<double>::AdditionalData());
-
-      InverseMatrix<SparseMatrix<double>,SparseILU<double> >
-      m_inverse (system_matrix.block(1,1), preconditioner);
-
-      // With the Schur complement and an
-      // efficient preconditioner at hand, we
-      // can solve the respective equation
-      // for the pressure (i.e. block 0 in
-      // the solution vector) in the usual
-      // way:
-      cg.solve (schur_complement, solution.block(1), schur_rhs,
-                m_inverse);
-
-      // After this first solution step, the
-      // hanging node constraints have to be
-      // distributed to the solution in order
-      // to achieve a consistent pressure
-      // field.
-      constraints.distribute (solution);
-
-      std::cout << "  "
-                << solver_control.last_step()
-                << " outer CG Schur complement iterations for pressure"
-                << std::endl;
-    }
-
-    // As in step-20, we finally need to
-    // solve for the velocity equation where
-    // we plug in the solution to the
-    // pressure equation. This involves only
-    // objects we already know - so we simply
-    // multiply $p$ by $B^T$, subtract the
-    // right hand side and multiply by the
-    // inverse of $A$. At the end, we need to
-    // distribute the constraints from
-    // hanging nodes in order to obtain a
-    // constistent flow field:
-    {
-      system_matrix.block(0,1).vmult (tmp, solution.block(1));
-      tmp *= -1;
-      tmp += system_rhs.block(0);
-
-      A_inverse.vmult (solution.block(0), tmp);
-
-      constraints.distribute (solution);
-    }
-  }
-
-
-  // @sect4{StokesProblem::output_results}
-
-  // The next function generates graphical
-  // output. In this example, we are going to
-  // use the VTK file format.  We attach
-  // names to the individual variables in the
-  // problem: <code>velocity</code> to the
-  // <code>dim</code> components of velocity
-  // and <code>pressure</code> to the
-  // pressure.
-  //
-  // Not all visualization programs have the
-  // ability to group individual vector
-  // components into a vector to provide
-  // vector plots; in particular, this holds
-  // for some VTK-based visualization
-  // programs. In this case, the logical
-  // grouping of components into vectors
-  // should already be described in the file
-  // containing the data. In other words,
-  // what we need to do is provide our output
-  // writers with a way to know which of the
-  // components of the finite element
-  // logically form a vector (with $d$
-  // components in $d$ space dimensions)
-  // rather than letting them assume that we
-  // simply have a bunch of scalar fields.
-  // This is achieved using the members of
-  // the
-  // <code>DataComponentInterpretation</code>
-  // namespace: as with the filename, we
-  // create a vector in which the first
-  // <code>dim</code> components refer to the
-  // velocities and are given the tag
-  // <code>DataComponentInterpretation::component_is_part_of_vector</code>;
-  // we finally push one tag
-  // <code>DataComponentInterpretation::component_is_scalar</code>
-  // to describe the grouping of the pressure
-  // variable.
-
-  // The rest of the function is then
-  // the same as in step-20.
-  template <int dim>
-  void
-  StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
-  {
-    std::vector<std::string> solution_names (dim, "velocity");
-    solution_names.push_back ("pressure");
-
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    data_component_interpretation
-    (dim, DataComponentInterpretation::component_is_part_of_vector);
-    data_component_interpretation
-    .push_back (DataComponentInterpretation::component_is_scalar);
-
-    DataOut<dim> data_out;
-    data_out.attach_dof_handler (dof_handler);
-    data_out.add_data_vector (solution, solution_names,
-                              DataOut<dim>::type_dof_data,
-                              data_component_interpretation);
-    data_out.build_patches ();
-
-    std::ostringstream filename;
-    filename << "solution-"
-             << Utilities::int_to_string (refinement_cycle, 2)
-             << ".vtk";
-
-    std::ofstream output (filename.str().c_str());
-    data_out.write_vtk (output);
-  }
-
-
-  // @sect4{StokesProblem::refine_mesh}
-
-  // This is the last interesting function of
-  // the <code>StokesProblem</code> class.
-  // As indicated by its name, it takes the
-  // solution to the problem and refines the
-  // mesh where this is needed. The procedure
-  // is the same as in the respective step in
-  // step-6, with the exception that we base
-  // the refinement only on the change in
-  // pressure, i.e., we call the Kelly error
-  // estimator with a mask
-  // object. Additionally, we do not coarsen
-  // the grid again:
-  template <int dim>
-  void
-  StokesProblem<dim>::refine_mesh ()
-  {
-    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-    std::vector<bool> component_mask (dim+1, false);
-    component_mask[dim] = true;
-    KellyErrorEstimator<dim>::estimate (dof_handler,
-                                        QGauss<dim-1>(degree+1),
-                                        typename FunctionMap<dim>::type(),
-                                        solution,
-                                        estimated_error_per_cell,
-                                        component_mask);
-
-    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                     estimated_error_per_cell,
-                                                     0.3, 0.0);
-    triangulation.execute_coarsening_and_refinement ();
-  }
-
-
-  // @sect4{StokesProblem::run}
-
-  // The last step in the Stokes class is, as
-  // usual, the function that generates the
-  // initial grid and calls the other
-  // functions in the respective order.
-  //
-  // We start off with a rectangle of size $4
-  // \times 1$ (in 2d) or $4 \times 1 \times
-  // 1$ (in 3d), placed in $R^2/R^3$ as
-  // $(-2,2)\times(-1,0)$ or
-  // $(-2,2)\times(0,1)\times(-1,0)$,
-  // respectively. It is natural to start
-  // with equal mesh size in each direction,
-  // so we subdivide the initial rectangle
-  // four times in the first coordinate
-  // direction. To limit the scope of the
-  // variables involved in the creation of
-  // the mesh to the range where we actually
-  // need them, we put the entire block
-  // between a pair of braces:
-  template <int dim>
-  void StokesProblem<dim>::run ()
-  {
-    {
-      std::vector<unsigned int> subdivisions (dim, 1);
-      subdivisions[0] = 4;
-
-      const Point<dim> bottom_left = (dim == 2 ?
-                                      Point<dim>(-2,-1) :
-                                      Point<dim>(-2,0,-1));
-      const Point<dim> top_right   = (dim == 2 ?
-                                      Point<dim>(2,0) :
-                                      Point<dim>(2,1,0));
-
-      GridGenerator::subdivided_hyper_rectangle (triangulation,
-                                                 subdivisions,
-                                                 bottom_left,
-                                                 top_right);
-    }
-
-    // A boundary indicator of 1 is set to all
-    // boundaries that are subject to Dirichlet
-    // boundary conditions, i.e.  to faces that
-    // are located at 0 in the last coordinate
-    // direction. See the example description
-    // above for details.
-    for (typename Triangulation<dim>::active_cell_iterator
-         cell = triangulation.begin_active();
-         cell != triangulation.end(); ++cell)
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->center()[dim-1] == 0)
-          cell->face(f)->set_all_boundary_ids(1);
-
-
-    // We then apply an initial refinement
-    // before solving for the first time. In
-    // 3D, there are going to be more degrees
-    // of freedom, so we refine less there:
-    triangulation.refine_global (4-dim);
-
-    // As first seen in step-6, we cycle over
-    // the different refinement levels and
-    // refine (except for the first cycle),
-    // setup the degrees of freedom and
-    // matrices, assemble, solve and create
-    // output:
-    for (unsigned int refinement_cycle = 0; refinement_cycle<6;
-         ++refinement_cycle)
-      {
-        std::cout << "Refinement cycle " << refinement_cycle << std::endl;
-
-        if (refinement_cycle > 0)
-          {
-            timer.enter_section("refine");
-            refine_mesh ();
-            timer.exit_section("refine");
-          }
-
-        timer.enter_section("setup");
-        setup_dofs ();
-        timer.exit_section("setup");
-
-        std::cout << "   Assembling..." << std::endl << std::flush;
-        timer.enter_section("assembly");
-        assemble_system ();
-        timer.exit_section("assembly");
-
-        std::cout << "   Solving..." << std::flush;
-        timer.enter_section("solver");
-        solve ();
-        timer.exit_section("solver");
-
-        timer.enter_section("results");
-        output_results (refinement_cycle);
-        timer.exit_section("results");
-
-        std::cout << std::endl;
-      }
-
-  }
-}
-
-
-// @sect3{The <code>main</code> function}
-
-// The main function is the same as in
-// step-20. We pass the element degree as a
-// parameter and choose the space dimension
-// at the well-known template slot.
-int main ()
-{
-  try
-    {
-      using namespace dealii;
-      using namespace Step22;
-
-      deallog.depth_console (0);
-
-      StokesProblem<2> flow_problem(1);
-      flow_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;
-
-      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/tests/benchmarks/sync.sh b/tests/benchmarks/sync.sh
deleted file mode 100755 (executable)
index 1bc58c4..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-#!/bin/bash
-echo "copying files to public_html/"
-cp *png ~/public_html/bench/
-cp index.html ~/public_html/bench/
-chmod a+r ~/public_html/bench/*png ~/public_html/bench/*html
-
diff --git a/tests/benchmarks/tablehandler/CMakeLists.txt b/tests/benchmarks/tablehandler/CMakeLists.txt
deleted file mode 100644 (file)
index 41f57c9..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-##
-#  CMake script for the step-1 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "table_handler")
-
-# Declare all source files the target consists of:
-SET(TARGET_SRC
-  ${TARGET}.cc
-  # You can specify additional files here!
-  )
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.0 QUIET
-  HINTS
-    ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR}
-  #
-  # If the deal.II library cannot be found (because it is not installed at a
-  # default location or your project resides at an uncommon place), you
-  # can specify additional hints for search paths here, e.g.
-  # "$ENV{HOME}/workspace/deal.II"
-  )
-
-IF (NOT ${deal.II_FOUND})
-   MESSAGE(FATAL_ERROR
-           "\n\n"
-          " *** Could not locate deal.II. *** "
-          "\n\n"
-           " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag 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/tests/benchmarks/tablehandler/table_handler.cc b/tests/benchmarks/tablehandler/table_handler.cc
deleted file mode 100644 (file)
index b5aaa66..0000000
+++ /dev/null
@@ -1,63 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2010 - 2014 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.
-//
-// ---------------------------------------------------------------------
-
-
-// we test the performance of writing a large table
-
-#include <deal.II/base/data_out_base.h>
-#include <deal.II/base/table_handler.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/timer.h>
-
-using namespace dealii;
-
-
-#include <vector>
-#include <iomanip>
-#include <iostream>
-#include <fstream>
-#include <string>
-
-
-int main ()
-{
-  TimerOutput timer (std::cout, TimerOutput::summary, TimerOutput::cpu_times);
-  deallog.depth_console(0);
-  deallog.threshold_double(1.e-10);
-
-  TableHandler table;
-
-  std::string keys[] = { "key1", "key2", "key3", "key4", "key5", "key6", "key7", "key8", "key9", "key10", "key11", "key12", "key13", "key14", "key15"};
-
-  unsigned int n_keys = 15;
-  unsigned int n_rows = 40000;
-
-  for (unsigned int j=0; j<n_rows; ++j)
-    {
-      table.add_value("begin", std::string("this is some text"));
-      for (unsigned int i=0; i<n_keys; ++i)
-        {
-          table.add_value(keys[i], j*1.0+i/100.0);
-        }
-    }
-
-  timer.enter_section("write");
-  {
-    std::ofstream data("datatable.txt");
-    table.write_text(data);
-//       TableHandler::table_with_separate_column_description);
-  }
-  timer.exit_section("write");
-}
diff --git a/tests/benchmarks/test_assembly/CMakeLists.txt b/tests/benchmarks/test_assembly/CMakeLists.txt
deleted file mode 100644 (file)
index f261836..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-##
-#  CMake script for the step-1 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "step-22")
-
-# Declare all source files the target consists of:
-SET(TARGET_SRC
-  ${TARGET}.cc
-  # You can specify additional files here!
-  )
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.0 QUIET
-  HINTS
-    ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR}
-  #
-  # If the deal.II library cannot be found (because it is not installed at a
-  # default location or your project resides at an uncommon place), you
-  # can specify additional hints for search paths here, e.g.
-  # "$ENV{HOME}/workspace/deal.II"
-  )
-
-IF (NOT ${deal.II_FOUND})
-   MESSAGE(FATAL_ERROR
-           "\n\n"
-          " *** Could not locate deal.II. *** "
-          "\n\n"
-           " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag 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/tests/benchmarks/test_assembly/step-22.cc b/tests/benchmarks/test_assembly/step-22.cc
deleted file mode 100644 (file)
index 33db9c9..0000000
+++ /dev/null
@@ -1,1339 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2008 - 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.
-//
-// ---------------------------------------------------------------------
-
-/*
- * step-22.cc
- */
-
-
-// @sect3{Include files}
-
-// As usual, we start by including
-// some well-known files:
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/constraint_matrix.h>
-
-#include <deal.II/grid/tria.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/tria_boundary_lib.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/grid_refinement.h>
-
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_renumbering.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_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 <deal.II/numerics/error_estimator.h>
-
-// Then we need to include the header file
-// for the sparse direct solver UMFPACK:
-#include <deal.II/lac/sparse_direct.h>
-#include <deal.II/base/timer.h>
-
-// This includes the library for the
-// incomplete LU factorization that will
-// be used as a preconditioner in 3D:
-#include <deal.II/lac/sparse_ilu.h>
-
-// This is C++:
-#include <fstream>
-#include <sstream>
-
-// As in all programs, the namespace dealii
-// is included:
-namespace Step22
-{
-  using namespace dealii;
-
-  // @sect3{Defining the inner preconditioner type}
-
-  // As explained in the introduction, we are
-  // going to use different preconditioners for
-  // two and three space dimensions,
-  // respectively. We distinguish between
-  // them by the use of the spatial dimension
-  // as a template parameter. See step-4 for
-  // details on templates. We are not going to
-  // create any preconditioner object here, all
-  // we do is to create class that holds a
-  // local typedef determining the
-  // preconditioner class so we can write our
-  // program in a dimension-independent way.
-  template <int dim>
-  struct InnerPreconditioner;
-
-  // In 2D, we are going to use a sparse direct
-  // solver as preconditioner:
-  template <>
-  struct InnerPreconditioner<2>
-  {
-    typedef SparseILU<double> type;
-//      typedef SparseDirectUMFPACK type;
-  };
-
-  // And the ILU preconditioning in 3D, called
-  // by SparseILU:
-  template <>
-  struct InnerPreconditioner<3>
-  {
-    typedef SparseILU<double> type;
-  };
-
-
-  // @sect3{The <code>StokesProblem</code> class template}
-
-  // This is an adaptation of step-20, so the
-  // main class and the data types are the
-  // same as used there. In this example we
-  // also use adaptive grid refinement, which
-  // is handled in analogy to
-  // step-6. According to the discussion in
-  // the introduction, we are also going to
-  // use the ConstraintMatrix for
-  // implementing Dirichlet boundary
-  // conditions. Hence, we change the name
-  // <code>hanging_node_constraints</code>
-  // into <code>constraints</code>.
-  template <int dim>
-  class StokesProblem
-  {
-  public:
-    StokesProblem (const unsigned int degree);
-    void run ();
-
-  private:
-    void setup_dofs ();
-    void assemble_system ();
-    void solve ();
-    void output_results (const unsigned int refinement_cycle) const;
-    void refine_mesh ();
-
-    const unsigned int   degree;
-
-    Triangulation<dim>   triangulation;
-    FESystem<dim>        fe;
-    DoFHandler<dim>      dof_handler;
-
-    ConstraintMatrix     constraints;
-
-    BlockSparsityPattern      sparsity_pattern;
-    BlockSparseMatrix<double> system_matrix;
-
-    BlockVector<double> solution;
-    BlockVector<double> system_rhs;
-
-    // This one is new: We shall use a
-    // so-called shared pointer structure to
-    // access the preconditioner. Shared
-    // pointers are essentially just a
-    // convenient form of pointers. Several
-    // shared pointers can point to the same
-    // object (just like regular pointers),
-    // but when the last shared pointer
-    // object to point to a preconditioner
-    // object is deleted (for example if a
-    // shared pointer object goes out of
-    // scope, if the class of which it is a
-    // member is destroyed, or if the pointer
-    // is assigned a different preconditioner
-    // object) then the preconditioner object
-    // pointed to is also destroyed. This
-    // ensures that we don't have to manually
-    // track in how many places a
-    // preconditioner object is still
-    // referenced, it can never create a
-    // memory leak, and can never produce a
-    // dangling pointer to an already
-    // destroyed object:
-    std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
-
-    TimerOutput timer;
-  };
-
-  // @sect3{Boundary values and right hand side}
-
-  // As in step-20 and most other
-  // example programs, the next task is
-  // to define the data for the PDE:
-  // For the Stokes problem, we are
-  // going to use natural boundary
-  // values on parts of the boundary
-  // (i.e. homogenous Neumann-type) for
-  // which we won't have to do anything
-  // special (the homogeneity implies
-  // that the corresponding terms in
-  // the weak form are simply zero),
-  // and boundary conditions on the
-  // velocity (Dirichlet-type) on the
-  // rest of the boundary, as described
-  // in the introduction.
-  //
-  // In order to enforce the Dirichlet
-  // boundary values on the velocity,
-  // we will use the
-  // VectorTools::interpolate_boundary_values
-  // function as usual which requires
-  // us to write a function object with
-  // as many components as the finite
-  // element has. In other words, we
-  // have to define the function on the
-  // $(u,p)$-space, but we are going to
-  // filter out the pressure component
-  // when interpolating the boundary
-  // values.
-
-  // The following function object is a
-  // representation of the boundary
-  // values described in the
-  // introduction:
-  template <int dim>
-  class BoundaryValues : public Function<dim>
-  {
-  public:
-    BoundaryValues () : Function<dim>(dim+1) {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
-  };
-
-
-  template <int dim>
-  double
-  BoundaryValues<dim>::value (const Point<dim>  &p,
-                              const unsigned int component) const
-  {
-    Assert (component < this->n_components,
-            ExcIndexRange (component, 0, this->n_components));
-
-    if (component == 0)
-      return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
-    return 0;
-  }
-
-
-  template <int dim>
-  void
-  BoundaryValues<dim>::vector_value (const Point<dim> &p,
-                                     Vector<double>   &values) const
-  {
-    for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = BoundaryValues<dim>::value (p, c);
-  }
-
-
-
-  // We implement similar functions for
-  // the right hand side which for the
-  // current example is simply zero:
-  template <int dim>
-  class RightHandSide : public Function<dim>
-  {
-  public:
-    RightHandSide () : Function<dim>(dim+1) {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
-
-  };
-
-
-  template <int dim>
-  double
-  RightHandSide<dim>::value (const Point<dim>  &/*p*/,
-                             const unsigned int /*component*/) const
-  {
-    return 0;
-  }
-
-
-  template <int dim>
-  void
-  RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                    Vector<double>   &values) const
-  {
-    for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = RightHandSide<dim>::value (p, c);
-  }
-
-
-  // @sect3{Linear solvers and preconditioners}
-
-  // The linear solvers and preconditioners are
-  // discussed extensively in the
-  // introduction. Here, we create the
-  // respective objects that will be used.
-
-  // @sect4{The <code>InverseMatrix</code> class template}
-
-  // The <code>InverseMatrix</code>
-  // class represents the data
-  // structure for an inverse
-  // matrix. It is derived from the one
-  // in step-20. The only difference is
-  // that we now do include a
-  // preconditioner to the matrix since
-  // we will apply this class to
-  // different kinds of matrices that
-  // will require different
-  // preconditioners (in step-20 we did
-  // not use a preconditioner in this
-  // class at all). The types of matrix
-  // and preconditioner are passed to
-  // this class via template
-  // parameters, and matrix and
-  // preconditioner objects of these
-  // types will then be passed to the
-  // constructor when an
-  // <code>InverseMatrix</code> object
-  // is created. The member function
-  // <code>vmult</code> is, as in
-  // step-20, a multiplication with a
-  // vector, obtained by solving a
-  // linear system:
-  template <class Matrix, class Preconditioner>
-  class InverseMatrix : public Subscriptor
-  {
-  public:
-    InverseMatrix (const Matrix         &m,
-                   const Preconditioner &preconditioner);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const Matrix> matrix;
-    const SmartPointer<const Preconditioner> preconditioner;
-  };
-
-
-  template <class Matrix, class Preconditioner>
-  InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
-                                                       const Preconditioner &preconditioner)
-    :
-    matrix (&m),
-    preconditioner (&preconditioner)
-  {}
-
-
-  // This is the implementation of the
-  // <code>vmult</code> function.
-
-  // In this class we use a rather large
-  // tolerance for the solver control. The
-  // reason for this is that the function is
-  // used very frequently, and hence, any
-  // additional effort to make the residual
-  // in the CG solve smaller makes the
-  // solution more expensive. Note that we do
-  // not only use this class as a
-  // preconditioner for the Schur complement,
-  // but also when forming the inverse of the
-  // Laplace matrix &ndash; which is hence
-  // directly responsible for the accuracy of
-  // the solution itself, so we can't choose
-  // a too large tolerance, either.
-  template <class Matrix, class Preconditioner>
-  void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
-                                                    const Vector<double> &src) const
-  {
-    SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
-    SolverCG<>    cg (solver_control);
-
-    dst = 0;
-
-    cg.solve (*matrix, dst, src, *preconditioner);
-  }
-
-
-  // @sect4{The <code>SchurComplement</code> class template}
-
-  // This class implements the Schur complement
-  // discussed in the introduction.  It is in
-  // analogy to step-20.  Though, we now call
-  // it with a template parameter
-  // <code>Preconditioner</code> in order to
-  // access that when specifying the respective
-  // type of the inverse matrix class. As a
-  // consequence of the definition above, the
-  // declaration <code>InverseMatrix</code> now
-  // contains the second template parameter
-  // for a preconditioner class as above, which
-  // affects the <code>SmartPointer</code>
-  // object <code>m_inverse</code> as well.
-  template <class Preconditioner>
-  class SchurComplement : public Subscriptor
-  {
-  public:
-    SchurComplement (const BlockSparseMatrix<double> &system_matrix,
-                     const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
-
-    mutable Vector<double> tmp1, tmp2;
-  };
-
-
-
-  template <class Preconditioner>
-  SchurComplement<Preconditioner>::
-  SchurComplement (const BlockSparseMatrix<double> &system_matrix,
-                   const InverseMatrix<SparseMatrix<double>,Preconditioner> &A_inverse)
-    :
-    system_matrix (&system_matrix),
-    A_inverse (&A_inverse),
-    tmp1 (system_matrix.block(0,0).m()),
-    tmp2 (system_matrix.block(0,0).m())
-  {}
-
-
-  template <class Preconditioner>
-  void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
-                                               const Vector<double> &src) const
-  {
-    system_matrix->block(0,1).vmult (tmp1, src);
-    A_inverse->vmult (tmp2, tmp1);
-    system_matrix->block(1,0).vmult (dst, tmp2);
-  }
-
-
-  // @sect3{StokesProblem class implementation}
-
-  // @sect4{StokesProblem::StokesProblem}
-
-  // The constructor of this class
-  // looks very similar to the one of
-  // step-20. The constructor
-  // initializes the variables for the
-  // polynomial degree, triangulation,
-  // finite element system and the dof
-  // handler. The underlying polynomial
-  // functions are of order
-  // <code>degree+1</code> for the
-  // vector-valued velocity components
-  // and of order <code>degree</code>
-  // for the pressure.  This gives the
-  // LBB-stable element pair
-  // $Q_{degree+1}^d\times Q_{degree}$,
-  // often referred to as the
-  // Taylor-Hood element.
-  //
-  // Note that we initialize the triangulation
-  // with a MeshSmoothing argument, which
-  // ensures that the refinement of cells is
-  // done in a way that the approximation of
-  // the PDE solution remains well-behaved
-  // (problems arise if grids are too
-  // unstructered), see the documentation of
-  // <code>Triangulation::MeshSmoothing</code>
-  // for details.
-  template <int dim>
-  StokesProblem<dim>::StokesProblem (const unsigned int degree)
-    :
-    degree (degree),
-    triangulation (Triangulation<dim>::maximum_smoothing),
-    fe (FE_Q<dim>(degree+1), dim,
-        FE_Q<dim>(degree), 1),
-    dof_handler (triangulation),
-    timer (std::cout, TimerOutput::summary, TimerOutput::cpu_times)
-  {}
-
-
-  // @sect4{StokesProblem::setup_dofs}
-
-  // Given a mesh, this function
-  // associates the degrees of freedom
-  // with it and creates the
-  // corresponding matrices and
-  // vectors. At the beginning it also
-  // releases the pointer to the
-  // preconditioner object (if the
-  // shared pointer pointed at anything
-  // at all at this point) since it
-  // will definitely not be needed any
-  // more after this point and will
-  // have to be re-computed after
-  // assembling the matrix, and unties
-  // the sparse matrix from its
-  // sparsity pattern object.
-  //
-  // We then proceed with distributing
-  // degrees of freedom and renumbering
-  // them: In order to make the ILU
-  // preconditioner (in 3D) work
-  // efficiently, it is important to
-  // enumerate the degrees of freedom
-  // in such a way that it reduces the
-  // bandwidth of the matrix, or maybe
-  // more importantly: in such a way
-  // that the ILU is as close as
-  // possible to a real LU
-  // decomposition. On the other hand,
-  // we need to preserve the block
-  // structure of velocity and pressure
-  // already seen in in step-20 and
-  // step-21. This is done in two
-  // steps: First, all dofs are
-  // renumbered to improve the ILU and
-  // then we renumber once again by
-  // components. Since
-  // <code>DoFRenumbering::component_wise</code>
-  // does not touch the renumbering
-  // within the individual blocks, the
-  // basic renumbering from the first
-  // step remains. As for how the
-  // renumber degrees of freedom to
-  // improve the ILU: deal.II has a
-  // number of algorithms that attempt
-  // to find orderings to improve ILUs,
-  // or reduce the bandwidth of
-  // matrices, or optimize some other
-  // aspect. The DoFRenumbering
-  // namespace shows a comparison of
-  // the results we obtain with several
-  // of these algorithms based on the
-  // testcase discussed here in this
-  // tutorial program. Here, we will
-  // use the traditional Cuthill-McKee
-  // algorithm already used in some of
-  // the previous tutorial programs.
-  // In the
-  // <a href="#improved-ilu">section on improved ILU</a>
-  // we're going to discuss this issue
-  // in more detail.
-
-  // There is one more change compared
-  // to previous tutorial programs:
-  // There is no reason in sorting the
-  // <code>dim</code> velocity
-  // components individually. In fact,
-  // rather than first enumerating all
-  // $x$-velocities, then all
-  // $y$-velocities, etc, we would like
-  // to keep all velocities at the same
-  // location together and only
-  // separate between velocities (all
-  // components) and pressures. By
-  // default, this is not what the
-  // DoFRenumbering::component_wise
-  // function does: it treats each
-  // vector component separately; what
-  // we have to do is group several
-  // components into "blocks" and pass
-  // this block structure to that
-  // function. Consequently, we
-  // allocate a vector
-  // <code>block_component</code> with
-  // as many elements as there are
-  // components and describe all
-  // velocity components to correspond
-  // to block 0, while the pressure
-  // component will form block 1:
-  template <int dim>
-  void StokesProblem<dim>::setup_dofs ()
-  {
-    A_preconditioner.reset ();
-    system_matrix.clear ();
-
-    timer.enter_section("distribute dofs");
-    dof_handler.distribute_dofs (fe);
-    timer.exit_section("distribute dofs");
-
-    timer.enter_section("renumbering");
-    DoFRenumbering::Cuthill_McKee (dof_handler);
-
-    std::vector<unsigned int> block_component (dim+1,0);
-    block_component[dim] = 1;
-    DoFRenumbering::component_wise (dof_handler, block_component);
-    timer.exit_section("renumbering");
-
-    timer.enter_section("setup constraints");
-    {
-      constraints.clear ();
-      std::vector<bool> component_mask (dim+1, true);
-      component_mask[dim] = false;
-      DoFTools::make_hanging_node_constraints (dof_handler,
-                                               constraints);
-      VectorTools::interpolate_boundary_values (dof_handler,
-                                                1,
-                                                BoundaryValues<dim>(),
-                                                constraints,
-                                                component_mask);
-    }
-
-    constraints.close ();
-    timer.exit_section("setup constraints");
-
-    // In analogy to step-20, we count the dofs
-    // in the individual components.  We could
-    // do this in the same way as there, but we
-    // want to operate on the block structure
-    // we used already for the renumbering: The
-    // function
-    // <code>DoFTools::count_dofs_per_block</code>
-    // does the same as
-    // <code>DoFTools::count_dofs_per_component</code>,
-    // but now grouped as velocity and pressure
-    // block via <code>block_component</code>.
-    std::vector<unsigned int> dofs_per_block (2);
-    DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
-    const unsigned int n_u = dofs_per_block[0],
-                       n_p = dofs_per_block[1];
-
-    std::cout << "   Number of active cells: "
-              << triangulation.n_active_cells()
-              << std::endl
-              << "   Number of degrees of freedom: "
-              << dof_handler.n_dofs()
-              << " (" << n_u << '+' << n_p << ')'
-              << std::endl;
-
-    // The next task is to allocate a
-    // sparsity pattern for the system matrix
-    // we will create. We could do this in
-    // the same way as in step-20,
-    // i.e. directly build an object of type
-    // SparsityPattern through
-    // DoFTools::make_sparsity_pattern. However,
-    // there is a major reason not to do so:
-    // In 3D, the function
-    // DoFTools::max_couplings_between_dofs
-    // yields a conservative but rather large
-    // number for the coupling between the
-    // individual dofs, so that the memory
-    // initially provided for the creation of
-    // the sparsity pattern of the matrix is
-    // far too much -- so much actually that
-    // the initial sparsity pattern won't
-    // even fit into the physical memory of
-    // most systems already for
-    // moderately-sized 3D problems, see also
-    // the discussion in step-18.  Instead,
-    // we first build a temporary object that
-    // uses a different data structure that
-    // doesn't require allocating more memory
-    // than necessary but isn't suitable for
-    // use as a basis of SparseMatrix or
-    // BlockSparseMatrix objects; in a second
-    // step we then copy this object into an
-    // object of BlockSparsityPattern. This
-    // is entirely analgous to what we
-    // already did in step-11 and step-18.
-    //
-    // The DynamicSparsityPattern class
-    // (and its block variant
-    // BlockDynamicSparsityPattern)
-    // has exactly the same interface, uses a
-    // different %internal data structure, and
-    // is linear in the number of degrees of
-    // freedom and therefore much more
-    // efficient for large problems.
-    //
-    // Consequently, we will use
-    // BlockDynamicSparsityPattern for our
-    // intermediate sparsity representation. All
-    // this is done inside a new scope, which
-    // means that the memory of <code>dsp</code>
-    // will be released once the information has
-    // been copied to
-    // <code>sparsity_pattern</code>.
-    {
-      timer.enter_section("make dsp");
-      BlockDynamicSparsityPattern dsp (2,2);
-
-      dsp.block(0,0).reinit (n_u, n_u);
-      dsp.block(1,0).reinit (n_p, n_u);
-      dsp.block(0,1).reinit (n_u, n_p);
-      dsp.block(1,1).reinit (n_p, n_p);
-
-      dsp.collect_sizes();
-
-      DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
-      timer.exit_section("make dsp");
-      timer.enter_section("copy sp");
-      sparsity_pattern.copy_from (dsp);
-      timer.exit_section("copy sp");
-    }
-
-    // Finally, the system matrix,
-    // solution and right hand side are
-    // created from the block
-    // structure as in step-20:
-    timer.enter_section("create matrix and vectors");
-    system_matrix.reinit (sparsity_pattern);
-
-    solution.reinit (2);
-    solution.block(0).reinit (n_u);
-    solution.block(1).reinit (n_p);
-    solution.collect_sizes ();
-
-    system_rhs.reinit (2);
-    system_rhs.block(0).reinit (n_u);
-    system_rhs.block(1).reinit (n_p);
-    system_rhs.collect_sizes ();
-    timer.exit_section("create matrix and vectors");
-  }
-
-
-  // @sect4{StokesProblem::assemble_system}
-
-  // The assembly process follows the
-  // discussion in step-20 and in the
-  // introduction. We use the well-known
-  // abbreviations for the data structures
-  // that hold the local matrix, right
-  // hand side, and global
-  // numbering of the degrees of freedom
-  // for the present cell.
-  template <int dim>
-  void StokesProblem<dim>::assemble_system ()
-  {
-    system_matrix=0;
-    system_rhs=0;
-
-    QGauss<dim>   quadrature_formula(degree+2);
-
-    FEValues<dim> fe_values (fe, quadrature_formula,
-                             update_values    |
-                             update_quadrature_points  |
-                             update_JxW_values |
-                             update_gradients);
-
-    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
-
-    const unsigned int   n_q_points      = quadrature_formula.size();
-
-    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>       local_rhs (dofs_per_cell);
-
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    const RightHandSide<dim>          right_hand_side;
-    std::vector<Vector<double> >      rhs_values (n_q_points,
-                                                  Vector<double>(dim+1));
-
-    // Next, we need two objects that work as
-    // extractors for the FEValues
-    // object. Their use is explained in detail
-    // in the report on @ref vector_valued :
-    const FEValuesExtractors::Vector velocities (0);
-    const FEValuesExtractors::Scalar pressure (dim);
-
-    // As an extension over step-20 and
-    // step-21, we include a few
-    // optimizations that make assembly
-    // much faster for this particular
-    // problem.  The improvements are
-    // based on the observation that we
-    // do a few calculations too many
-    // times when we do as in step-20:
-    // The symmetric gradient actually
-    // has <code>dofs_per_cell</code>
-    // different values per quadrature
-    // point, but we extract it
-    // <code>dofs_per_cell*dofs_per_cell</code>
-    // times from the FEValues object -
-    // for both the loop over
-    // <code>i</code> and the inner
-    // loop over <code>j</code>. In 3d,
-    // that means evaluating it
-    // $89^2=7921$ instead of $89$
-    // times, a not insignificant
-    // difference.
-    //
-    // So what we're
-    // going to do here is to avoid
-    // such repeated calculations by
-    // getting a vector of rank-2
-    // tensors (and similarly for
-    // the divergence and the basis
-    // function value on pressure)
-    // at the quadrature point prior
-    // to starting the loop over the
-    // dofs on the cell. First, we
-    // create the respective objects
-    // that will hold these
-    // values. Then, we start the
-    // loop over all cells and the loop
-    // over the quadrature points,
-    // where we first extract these
-    // values. There is one more
-    // optimization we implement here:
-    // the local matrix (as well as
-    // the global one) is going to
-    // be symmetric, since all
-    // the operations involved are
-    // symmetric with respect to $i$
-    // and $j$. This is implemented by
-    // simply running the inner loop
-    // not to <code>dofs_per_cell</code>,
-    // but only up to <code>i</code>,
-    // the index of the outer loop.
-    std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
-    std::vector<double>                  div_phi_u   (dofs_per_cell);
-    std::vector<double>                  phi_p       (dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-    for (; cell!=endc; ++cell)
-      {
-        fe_values.reinit (cell);
-        local_matrix = 0;
-        local_rhs = 0;
-
-        right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
-                                          rhs_values);
-
-        for (unsigned int q=0; q<n_q_points; ++q)
-          {
-            for (unsigned int k=0; k<dofs_per_cell; ++k)
-              {
-                symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
-                div_phi_u[k]     = fe_values[velocities].divergence (k, q);
-                phi_p[k]         = fe_values[pressure].value (k, q);
-              }
-
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              {
-                for (unsigned int j=0; j<=i; ++j)
-                  {
-                    local_matrix(i,j) += (2 * (symgrad_phi_u[i] * symgrad_phi_u[j])
-                                          - div_phi_u[i] * phi_p[j]
-                                          - phi_p[i] * div_phi_u[j]
-                                          + phi_p[i] * phi_p[j])
-                                         * fe_values.JxW(q);
-
-                  }
-
-                const unsigned int component_i =
-                  fe.system_to_component_index(i).first;
-                local_rhs(i) += fe_values.shape_value(i,q) *
-                                rhs_values[q](component_i) *
-                                fe_values.JxW(q);
-              }
-          }
-
-        // Note that in the above computation
-        // of the local matrix contribution
-        // we added the term <code> phi_p[i] *
-        // phi_p[j] </code>, yielding a
-        // pressure mass matrix in the
-        // $(1,1)$ block of the matrix as
-        // discussed in the
-        // introduction. That this term only
-        // ends up in the $(1,1)$ block stems
-        // from the fact that both of the
-        // factors in <code>phi_p[i] *
-        // phi_p[j]</code> are only non-zero
-        // when all the other terms vanish
-        // (and the other way around).
-        //
-        // Note also that operator* is
-        // overloaded for symmetric
-        // tensors, yielding the scalar
-        // product between the two
-        // tensors in the first line of
-        // the local matrix
-        // contribution.
-
-        // Before we can write the local data
-        // into the global matrix (and
-        // simultaneously use the
-        // ConstraintMatrix object to apply
-        // Dirichlet boundary conditions and
-        // eliminate hanging node
-        // constraints, as we discussed in
-        // the introduction), we have to be
-        // careful about one thing,
-        // though. We have only build up half
-        // of the local matrix because of
-        // symmetry, but we're going to save
-        // the full system matrix in order to
-        // use the standard functions for
-        // solution. This is done by flipping
-        // the indices in case we are
-        // pointing into the empty part of
-        // the local matrix.
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          for (unsigned int j=i+1; j<dofs_per_cell; ++j)
-            local_matrix(i,j) = local_matrix(j,i);
-
-        cell->get_dof_indices (local_dof_indices);
-        constraints.distribute_local_to_global (local_matrix, local_rhs,
-                                                local_dof_indices,
-                                                system_matrix, system_rhs);
-      }
-
-    // Before we're going to solve this
-    // linear system, we generate a
-    // preconditioner for the
-    // velocity-velocity matrix, i.e.,
-    // <code>block(0,0)</code> in the
-    // system matrix. As mentioned
-    // above, this depends on the
-    // spatial dimension. Since the two
-    // classes described by the
-    // <code>InnerPreconditioner::type</code>
-    // typedef have the same interface,
-    // we do not have to do anything
-    // different whether we want to use
-    // a sparse direct solver or an
-    // ILU:
-    std::cout << "   Computing preconditioner..." << std::endl << std::flush;
-
-    A_preconditioner
-      = std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
-    A_preconditioner->initialize (system_matrix.block(0,0),
-                                  typename InnerPreconditioner<dim>::type::AdditionalData());
-
-  }
-
-
-
-  // @sect4{StokesProblem::solve}
-
-  // After the discussion in the introduction
-  // and the definition of the respective
-  // classes above, the implementation of the
-  // <code>solve</code> function is rather
-  // straigt-forward and done in a similar way
-  // as in step-20. To start with, we need an
-  // object of the <code>InverseMatrix</code>
-  // class that represents the inverse of the
-  // matrix A. As described in the
-  // introduction, the inverse is generated
-  // with the help of an inner preconditioner
-  // of type
-  // <code>InnerPreconditioner::type</code>.
-  template <int dim>
-  void StokesProblem<dim>::solve ()
-  {
-    const InverseMatrix<SparseMatrix<double>,
-          typename InnerPreconditioner<dim>::type>
-          A_inverse (system_matrix.block(0,0), *A_preconditioner);
-    Vector<double> tmp (solution.block(0).size());
-
-    // This is as in step-20. We generate the
-    // right hand side $B A^{-1} F - G$ for the
-    // Schur complement and an object that
-    // represents the respective linear
-    // operation $B A^{-1} B^T$, now with a
-    // template parameter indicating the
-    // preconditioner - in accordance with the
-    // definition of the class.
-    {
-      Vector<double> schur_rhs (solution.block(1).size());
-      A_inverse.vmult (tmp, system_rhs.block(0));
-      system_matrix.block(1,0).vmult (schur_rhs, tmp);
-      schur_rhs -= system_rhs.block(1);
-
-      SchurComplement<typename InnerPreconditioner<dim>::type>
-      schur_complement (system_matrix, A_inverse);
-
-      // The usual control structures for
-      // the solver call are created...
-      SolverControl solver_control (solution.block(1).size(),
-                                    1e-6*schur_rhs.l2_norm());
-      SolverCG<>    cg (solver_control);
-
-      // Now to the preconditioner to the
-      // Schur complement. As explained in
-      // the introduction, the
-      // preconditioning is done by a mass
-      // matrix in the pressure variable.  It
-      // is stored in the $(1,1)$ block of
-      // the system matrix (that is not used
-      // anywhere else but in
-      // preconditioning).
-      //
-      // Actually, the solver needs to have
-      // the preconditioner in the form
-      // $P^{-1}$, so we need to create an
-      // inverse operation. Once again, we
-      // use an object of the class
-      // <code>InverseMatrix</code>, which
-      // implements the <code>vmult</code>
-      // operation that is needed by the
-      // solver.  In this case, we have to
-      // invert the pressure mass matrix. As
-      // it already turned out in earlier
-      // tutorial programs, the inversion of
-      // a mass matrix is a rather cheap and
-      // straight-forward operation (compared
-      // to, e.g., a Laplace matrix). The CG
-      // method with ILU preconditioning
-      // converges in 5-10 steps,
-      // independently on the mesh size.
-      // This is precisely what we do here:
-      // We choose another ILU preconditioner
-      // and take it along to the
-      // InverseMatrix object via the
-      // corresponding template parameter.  A
-      // CG solver is then called within the
-      // vmult operation of the inverse
-      // matrix.
-      //
-      // An alternative that is cheaper to
-      // build, but needs more iterations
-      // afterwards, would be to choose a
-      // SSOR preconditioner with factor
-      // 1.2. It needs about twice the number
-      // of iterations, but the costs for its
-      // generation are almost neglible.
-      SparseILU<double> preconditioner;
-      preconditioner.initialize (system_matrix.block(1,1),
-                                 SparseILU<double>::AdditionalData());
-
-      InverseMatrix<SparseMatrix<double>,SparseILU<double> >
-      m_inverse (system_matrix.block(1,1), preconditioner);
-
-      // With the Schur complement and an
-      // efficient preconditioner at hand, we
-      // can solve the respective equation
-      // for the pressure (i.e. block 0 in
-      // the solution vector) in the usual
-      // way:
-      cg.solve (schur_complement, solution.block(1), schur_rhs,
-                m_inverse);
-
-      // After this first solution step, the
-      // hanging node constraints have to be
-      // distributed to the solution in order
-      // to achieve a consistent pressure
-      // field.
-      constraints.distribute (solution);
-
-      std::cout << "  "
-                << solver_control.last_step()
-                << " outer CG Schur complement iterations for pressure"
-                << std::endl;
-    }
-
-    // As in step-20, we finally need to
-    // solve for the velocity equation where
-    // we plug in the solution to the
-    // pressure equation. This involves only
-    // objects we already know - so we simply
-    // multiply $p$ by $B^T$, subtract the
-    // right hand side and multiply by the
-    // inverse of $A$. At the end, we need to
-    // distribute the constraints from
-    // hanging nodes in order to obtain a
-    // constistent flow field:
-    {
-      system_matrix.block(0,1).vmult (tmp, solution.block(1));
-      tmp *= -1;
-      tmp += system_rhs.block(0);
-
-      A_inverse.vmult (solution.block(0), tmp);
-
-      constraints.distribute (solution);
-    }
-  }
-
-
-  // @sect4{StokesProblem::output_results}
-
-  // The next function generates graphical
-  // output. In this example, we are going to
-  // use the VTK file format.  We attach
-  // names to the individual variables in the
-  // problem: <code>velocity</code> to the
-  // <code>dim</code> components of velocity
-  // and <code>pressure</code> to the
-  // pressure.
-  //
-  // Not all visualization programs have the
-  // ability to group individual vector
-  // components into a vector to provide
-  // vector plots; in particular, this holds
-  // for some VTK-based visualization
-  // programs. In this case, the logical
-  // grouping of components into vectors
-  // should already be described in the file
-  // containing the data. In other words,
-  // what we need to do is provide our output
-  // writers with a way to know which of the
-  // components of the finite element
-  // logically form a vector (with $d$
-  // components in $d$ space dimensions)
-  // rather than letting them assume that we
-  // simply have a bunch of scalar fields.
-  // This is achieved using the members of
-  // the
-  // <code>DataComponentInterpretation</code>
-  // namespace: as with the filename, we
-  // create a vector in which the first
-  // <code>dim</code> components refer to the
-  // velocities and are given the tag
-  // <code>DataComponentInterpretation::component_is_part_of_vector</code>;
-  // we finally push one tag
-  // <code>DataComponentInterpretation::component_is_scalar</code>
-  // to describe the grouping of the pressure
-  // variable.
-
-  // The rest of the function is then
-  // the same as in step-20.
-  template <int dim>
-  void
-  StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
-  {
-    std::vector<std::string> solution_names (dim, "velocity");
-    solution_names.push_back ("pressure");
-
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    data_component_interpretation
-    (dim, DataComponentInterpretation::component_is_part_of_vector);
-    data_component_interpretation
-    .push_back (DataComponentInterpretation::component_is_scalar);
-
-    DataOut<dim> data_out;
-    data_out.attach_dof_handler (dof_handler);
-    data_out.add_data_vector (solution, solution_names,
-                              DataOut<dim>::type_dof_data,
-                              data_component_interpretation);
-    data_out.build_patches ();
-
-    std::ostringstream filename;
-    filename << "solution-"
-             << Utilities::int_to_string (refinement_cycle, 2)
-             << ".vtk";
-
-    std::ofstream output (filename.str().c_str());
-    data_out.write_vtk (output);
-  }
-
-
-  // @sect4{StokesProblem::refine_mesh}
-
-  // This is the last interesting function of
-  // the <code>StokesProblem</code> class.
-  // As indicated by its name, it takes the
-  // solution to the problem and refines the
-  // mesh where this is needed. The procedure
-  // is the same as in the respective step in
-  // step-6, with the exception that we base
-  // the refinement only on the change in
-  // pressure, i.e., we call the Kelly error
-  // estimator with a mask
-  // object. Additionally, we do not coarsen
-  // the grid again:
-  template <int dim>
-  void
-  StokesProblem<dim>::refine_mesh ()
-  {
-    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-    srand (0);
-    for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
-      estimated_error_per_cell(i) = 1.0*(myrand() % 1000);
-    /*
-        std::vector<bool> component_mask (dim+1, false);
-        component_mask[dim] = true;
-        KellyErrorEstimator<dim>::estimate (dof_handler,
-                                            QGauss<dim-1>(degree+1),
-                                            typename FunctionMap<dim>::type(),
-                                            solution,
-                                            estimated_error_per_cell,
-                                            component_mask);
-    */
-    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                     estimated_error_per_cell,
-                                                     0.3, 0.0);
-    triangulation.execute_coarsening_and_refinement ();
-  }
-
-
-  // @sect4{StokesProblem::run}
-
-  // The last step in the Stokes class is, as
-  // usual, the function that generates the
-  // initial grid and calls the other
-  // functions in the respective order.
-  //
-  // We start off with a rectangle of size $4
-  // \times 1$ (in 2d) or $4 \times 1 \times
-  // 1$ (in 3d), placed in $R^2/R^3$ as
-  // $(-2,2)\times(-1,0)$ or
-  // $(-2,2)\times(0,1)\times(-1,0)$,
-  // respectively. It is natural to start
-  // with equal mesh size in each direction,
-  // so we subdivide the initial rectangle
-  // four times in the first coordinate
-  // direction. To limit the scope of the
-  // variables involved in the creation of
-  // the mesh to the range where we actually
-  // need them, we put the entire block
-  // between a pair of braces:
-  template <int dim>
-  void StokesProblem<dim>::run ()
-  {
-    {
-      std::vector<unsigned int> subdivisions (dim, 1);
-      subdivisions[0] = 4;
-
-      const Point<dim> bottom_left = (dim == 2 ?
-                                      Point<dim>(-2,-1) :
-                                      Point<dim>(-2,0,-1));
-      const Point<dim> top_right   = (dim == 2 ?
-                                      Point<dim>(2,0) :
-                                      Point<dim>(2,1,0));
-
-      GridGenerator::subdivided_hyper_rectangle (triangulation,
-                                                 subdivisions,
-                                                 bottom_left,
-                                                 top_right);
-    }
-
-    // A boundary indicator of 1 is set to all
-    // boundaries that are subject to Dirichlet
-    // boundary conditions, i.e.  to faces that
-    // are located at 0 in the last coordinate
-    // direction. See the example description
-    // above for details.
-    for (typename Triangulation<dim>::active_cell_iterator
-         cell = triangulation.begin_active();
-         cell != triangulation.end(); ++cell)
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->center()[dim-1] == 0)
-          cell->face(f)->set_all_boundary_ids(1);
-
-
-    // We then apply an initial refinement
-    // before solving for the first time. In
-    // 3D, there are going to be more degrees
-    // of freedom, so we refine less there:
-    triangulation.refine_global (8-dim);
-
-    // As first seen in step-6, we cycle over
-    // the different refinement levels and
-    // refine (except for the first cycle),
-    // setup the degrees of freedom and
-    // matrices, assemble, solve and create
-    // output:
-    for (unsigned int refinement_cycle = 0; refinement_cycle<1;
-         ++refinement_cycle)
-      {
-        std::cout << "Refinement cycle " << refinement_cycle << std::endl;
-
-        //if (refinement_cycle > 0)
-        {
-          timer.enter_section("refine");
-          refine_mesh ();
-          refine_mesh ();
-          refine_mesh ();
-          //triangulation.refine_global(1);
-          timer.exit_section("refine");
-        }
-
-        //timer.enter_section("setup");
-        setup_dofs ();
-        //timer.exit_section("setup");
-
-        std::cout << "   Assembling..." << std::endl << std::flush;
-        timer.enter_section("assembly");
-        assemble_system ();
-        timer.exit_section("assembly");
-
-        std::cout << "   Solving..." << std::flush;
-        //timer.enter_section("solver");
-        //solve ();
-        //timer.exit_section("solver");
-
-        //timer.enter_section("results");
-        //output_results (refinement_cycle);
-        //timer.exit_section("results");
-
-        std::cout << std::endl;
-      }
-
-  }
-}
-
-
-// @sect3{The <code>main</code> function}
-
-// The main function is the same as in
-// step-20. We pass the element degree as a
-// parameter and choose the space dimension
-// at the well-known template slot.
-int main ()
-{
-  try
-    {
-      using namespace dealii;
-      using namespace Step22;
-
-      deallog.depth_console (0);
-
-      StokesProblem<2> flow_problem(1);
-      flow_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;
-
-      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/tests/benchmarks/test_hp/CMakeLists.txt b/tests/benchmarks/test_hp/CMakeLists.txt
deleted file mode 100644 (file)
index 0d9a9b4..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-##
-#  CMake script for the step-1 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "step")
-
-# Declare all source files the target consists of:
-SET(TARGET_SRC
-  ${TARGET}.cc
-  # You can specify additional files here!
-  )
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.0 QUIET
-  HINTS
-    ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR}
-  #
-  # If the deal.II library cannot be found (because it is not installed at a
-  # default location or your project resides at an uncommon place), you
-  # can specify additional hints for search paths here, e.g.
-  # "$ENV{HOME}/workspace/deal.II"
-  )
-
-IF (NOT ${deal.II_FOUND})
-   MESSAGE(FATAL_ERROR
-           "\n\n"
-          " *** Could not locate deal.II. *** "
-          "\n\n"
-           " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag 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/tests/benchmarks/test_hp/step.cc b/tests/benchmarks/test_hp/step.cc
deleted file mode 100644 (file)
index 0e47d46..0000000
+++ /dev/null
@@ -1,931 +0,0 @@
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2000 - 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
- */
-
-// @sect3{Include files}
-
-// The first few files have already been covered in previous examples and will
-// thus not be further commented on.
-#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/compressed_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/tria_boundary_lib.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.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 <deal.II/hp/dof_handler.h>
-#include <deal.II/hp/fe_values.h>
-#include <deal.II/base/timer.h>
-
-#include <deal.II/lac/sparse_direct.h>
-
-#include <fstream>
-#include <iostream>
-
-// From the following include file we will import the declaration of
-// H1-conforming finite element shape functions. This family of finite
-// elements is called <code>FE_Q</code>, and was used in all examples before
-// already to define the usual bi- or tri-linear elements, but we will now use
-// it for bi-quadratic elements:
-#include <deal.II/fe/fe_q.h>
-// We will not read the grid from a file as in the previous example, but
-// generate it using a function of the library. However, we will want to write
-// out the locally refined grids (just the grid, not the solution) in each
-// step, so we need the following include file instead of
-// <code>grid_in.h</code>:
-#include <deal.II/grid/grid_out.h>
-
-
-// When using locally refined grids, we will get so-called <code>hanging
-// nodes</code>. However, the standard finite element methods assumes that the
-// discrete solution spaces be continuous, so we need to make sure that the
-// degrees of freedom on hanging nodes conform to some constraints such that
-// the global solution is continuous. We are also going to store the boundary
-// conditions in this object. The following file contains a class which is
-// used to handle these constraints:
-#include <deal.II/lac/constraint_matrix.h>
-
-// In order to refine our grids locally, we need a function from the library
-// that decides which cells to flag for refinement or coarsening based on the
-// error indicators we have computed. This function is defined here:
-#include <deal.II/grid/grid_refinement.h>
-
-// Finally, we need a simple way to actually compute the refinement indicators
-// based on some error estimat. While in general, adaptivity is very
-// problem-specific, the error indicator in the following file often yields
-// quite nicely adapted grids for a wide class of problems.
-#include <deal.II/numerics/error_estimator.h>
-
-// Finally, this is as in previous programs:
-using namespace dealii;
-
-
-// @sect3{The <code>Step6</code> class template}
-
-// The main class is again almost unchanged. Two additions, however, are made:
-// we have added the <code>refine_grid</code> function, which is used to
-// adaptively refine the grid (instead of the global refinement in the
-// previous examples), and a variable which will hold the constraints. In
-// addition, we have added a destructor to the class for reasons that will
-// become clear when we discuss its implementation.
-
-template <int dim>
-class Step6
-{
-public:
-  Step6 ();
-  ~Step6 ();
-
-  void run ();
-
-private:
-  void setup_system ();
-  void setup_system_hp ();
-  void assemble_system ();
-  void assemble_system_hp ();
-  void solve ();
-  void refine_grid ();
-  void output_results (const unsigned int cycle) const;
-
-  Triangulation<dim>   triangulation;
-
-  DoFHandler<dim>      dof_handler;
-  hp::DoFHandler<dim>      hpdof_handler;
-  FE_Q<dim> fe;
-  hp::FECollection<dim> hpfe;
-
-  // This is the new variable in the main class. We need an object which holds
-  // a list of constraints to hold the hanging nodes and the boundary
-  // conditions.
-  ConstraintMatrix     constraints;
-  TimerOutput                               computing_timer;
-
-  SparsityPattern      sparsity_pattern;
-  SparseMatrix<double> system_matrix;
-
-  Vector<double>       solution;
-  Vector<double>       system_rhs;
-};
-
-
-
-// @sect3{Nonconstant coefficients}
-
-// The implementation of nonconstant coefficients is copied verbatim from
-// step-5:
-
-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;
-    }
-}
-
-
-// @sect3{The <code>Step6</code> class implementation}
-
-// @sect4{Step6::Step6}
-
-// The constructor of this class is mostly the same as before, but this time
-// we want to use the quadratic element. To do so, we only have to replace the
-// constructor argument (which was <code>1</code> in all previous examples) by
-// the desired polynomial degree (here <code>2</code>):
-template <int dim>
-Step6<dim>::Step6 ()
-  :
-  dof_handler (triangulation),
-  hpdof_handler (triangulation),
-  fe (3),
-  computing_timer (std::cout,
-                   TimerOutput::summary,
-                   TimerOutput::wall_times)
-
-{
-  hpfe.push_back (FE_Q<dim>(3));
-}
-
-
-// @sect4{Step6::~Step6}
-
-// Here comes the added destructor of the class. The reason why we want to add
-// it is a subtle change in the order of data elements in the class as
-// compared to all previous examples: the <code>dof_handler</code> object was
-// defined before and not after the <code>fe</code> object. Of course we could
-// have left this order unchanged, but we would like to show what happens if
-// the order is reversed since this produces a rather nasty side-effect and
-// results in an error which is difficult to track down if one does not know
-// what happens.
-//
-// Basically what happens is the following: when we distribute the degrees of
-// freedom using the function call <code>dof_handler.distribute_dofs()</code>,
-// the <code>dof_handler</code> also stores a pointer to the finite element in
-// use. Since this pointer is used every now and then until either the degrees
-// of freedom are re-distributed using another finite element object or until
-// the <code>dof_handler</code> object is destroyed, it would be unwise if we
-// would allow the finite element object to be deleted before the
-// <code>dof_handler</code> object. To disallow this, the DoF handler
-// increases a counter inside the finite element object which counts how many
-// objects use that finite element (this is what the
-// <code>Subscriptor</code>/<code>SmartPointer</code> class pair is used for,
-// in case you want something like this for your own programs; see step-7 for
-// a more complete discussion of this topic). The finite element object will
-// refuse its destruction if that counter is larger than zero, since then some
-// other objects might rely on the persistence of the finite element
-// object. An exception will then be thrown and the program will usually abort
-// upon the attempt to destroy the finite element.
-//
-// To be fair, such exceptions about still used objects are not particularly
-// popular among programmers using deal.II, since they only tell us that
-// something is wrong, namely that some other object is still using the object
-// that is presently being destructed, but most of the time not who this user
-// is. It is therefore often rather time-consuming to find out where the
-// problem exactly is, although it is then usually straightforward to remedy
-// the situation. However, we believe that the effort to find invalid
-// references to objects that do no longer exist is less if the problem is
-// detected once the reference becomes invalid, rather than when non-existent
-// objects are actually accessed again, since then usually only invalid data
-// is accessed, but no error is immediately raised.
-//
-// Coming back to the present situation, if we did not write this destructor,
-// the compiler will generate code that triggers exactly the behavior sketched
-// above. The reason is that member variables of the <code>Step6</code> class
-// are destructed bottom-up (i.e. in reverse order of their declaration in the
-// class), as always in C++. Thus, the finite element object will be
-// destructed before the DoF handler object, since its declaration is below
-// the one of the DoF handler. This triggers the situation above, and an
-// exception will be raised when the <code>fe</code> object is
-// destructed. What needs to be done is to tell the <code>dof_handler</code>
-// object to release its lock to the finite element. Of course, the
-// <code>dof_handler</code> will only release its lock if it really does not
-// need the finite element any more, i.e. when all finite element related data
-// is deleted from it. For this purpose, the <code>DoFHandler</code> class has
-// a function <code>clear</code> which deletes all degrees of freedom, and
-// releases its lock to the finite element. After this, you can safely
-// destruct the finite element object since its internal counter is then zero.
-//
-// For completeness, we add the output of the exception that would have been
-// triggered without this destructor, to the end of the results section of
-// this example.
-template <int dim>
-Step6<dim>::~Step6 ()
-{
-  dof_handler.clear ();
-}
-
-
-// @sect4{Step6::setup_system}
-
-// The next function is setting up all the variables that describe the linear
-// finite element problem, such as the DoF handler, the matrices, and
-// vectors. The difference to what we did in step-5 is only that we now also
-// have to take care of handing node constraints. These constraints are
-// handled almost transparently by the library, i.e. you only need to know
-// that they exist and how to get them, but you do not have to know how they
-// are formed or what exactly is done with them.
-//
-// At the beginning of the function, you find all the things that are the same
-// as in step-5: setting up the degrees of freedom (this time we have
-// quadratic elements, but there is no difference from a user code perspective
-// to the linear -- or cubic, for that matter -- case), generating the
-// sparsity pattern, and initializing the solution and right hand side
-// vectors. Note that the sparsity pattern will have significantly more
-// entries per row now, since there are now 9 degrees of freedom per cell, not
-// only four, that can couple with each other. The
-// <code>dof_Handler.max_couplings_between_dofs()</code> call will take care
-// of this, however:
-template <int dim>
-void Step6<dim>::setup_system ()
-{
-  computing_timer.enter_section ("distribute");
-  dof_handler.distribute_dofs (fe);
-  computing_timer.exit_section ("distribute");
-
-  computing_timer.enter_section ("distribute_hp");
-  hpdof_handler.distribute_dofs (hpfe);
-  computing_timer.exit_section ("distribute_hp");
-
-  solution.reinit (dof_handler.n_dofs());
-  system_rhs.reinit (dof_handler.n_dofs());
-
-
-  // After setting up all the degrees of freedoms, here are now the
-  // differences compared to step-5, all of which are related to constraints
-  // associated with the hanging nodes. In the class desclaration, we have
-  // already allocated space for an object <code>constraints</code> that will
-  // hold a list of these constraints (they form a matrix, which is reflected
-  // in the name of the class, but that is immaterial for the moment). Now we
-  // have to fill this object. This is done using the following function calls
-  // (the first clears the contents of the object that may still be left over
-  // from computations on the previous mesh before the last adaptive
-  // refinement):
-  constraints.clear ();
-  //computing_timer.enter_section ("hanging");
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                           constraints);
-  //computing_timer.exit_section ("hanging");
-  /*  constraints.clear ();
-      computing_timer.enter_section ("hanging_hp");
-    DoFTools::make_hanging_node_constraints (dof_handler,
-                                             constraints);
-      computing_timer.exit_section ("hanging_hp");
-  */
-
-  // Now we are ready to interpolate the ZeroFunction to our boundary with
-  // indicator 0 (the whole boundary) and store the resulting constraints in
-  // our <code>constraints</code> object. Note that we do not to apply the
-  // boundary conditions after assembly, like we did in earlier steps.  As
-  // almost all the stuff, the interpolation of boundary values works also for
-  // higher order elements without the need to change your code for that. We
-  // note that for proper results, it is important that the elimination of
-  // boundary nodes from the system of equations happens *after* the
-  // elimination of hanging nodes. For that reason we are filling the boundary
-  // values into the ContraintMatrix after the hanging node constraints.
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                            0,
-                                            ZeroFunction<dim>(),
-                                            constraints);
-
-
-  // The next step is <code>closing</code> this object. After all constraints
-  // have been added, they need to be sorted and rearranged to perform some
-  // actions more efficiently. This postprocessing is done using the
-  // <code>close()</code> function, after which no further constraints may be
-  // added any more:
-  constraints.close ();
-
-  // Now we first build our compressed sparsity pattern like we did in the
-  // previous examples. Nevertheless, we do not copy it to the final sparsity
-  // pattern immediately.  Note that we call a variant of
-  // make_sparsity_pattern that takes the ConstraintMatrix as the third
-  // argument. We are letting the routine know that we will never write into
-  // the locations given by <code>constraints</code> by setting the argument
-  // <code>keep_constrained_dofs</code> to false (in other words, that we will
-  // never write into entries of the matrix that correspond to constrained
-  // degrees of freedom). If we were to condense the
-  // constraints after assembling, we would have to pass <code>true</code>
-  // instead because then we would first write into these locations only to
-  // later set them to zero again during condensation.
-  CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
-//    computing_timer.enter_section ("makesp");
-  DoFTools::make_sparsity_pattern(dof_handler,
-                                  c_sparsity,
-                                  constraints,
-                                  /*keep_constrained_dofs = */ false);
-//    computing_timer.exit_section ("makesp");
-
-  // Now all non-zero entries of the matrix are known (i.e. those from
-  // regularly assembling the matrix and those that were introduced by
-  // eliminating constraints). We can thus copy our intermediate object to the
-  // sparsity pattern:
-  sparsity_pattern.copy_from(c_sparsity);
-
-  // Finally, the so-constructed sparsity pattern serves as the basis on top
-  // of which we will create the sparse matrix:
-  system_matrix.reinit (sparsity_pattern);
-}
-
-// @sect4{Step6::assemble_system}
-
-// Next, we have to assemble the matrix again. There are two code changes
-// compared to step-5:
-//
-// First, we have to use a higher-order quadrature formula to account for the
-// higher polynomial degree in the finite element shape functions. This is
-// easy to change: the constructor of the <code>QGauss</code> class takes the
-// number of quadrature points in each space direction. Previously, we had two
-// points for bilinear elements. Now we should use three points for
-// biquadratic elements.
-//
-// Second, to copy the local matrix and vector on each cell into the global
-// system, we are no longer using a hand-written loop. Instead, we use
-// <code>ConstraintMatrix::distribute_local_to_global</code> that internally
-// executes this loop and eliminates all the constraints at the same time.
-//
-// The rest of the code that forms the local contributions remains
-// unchanged. It is worth noting, however, that under the hood several things
-// are different than before. First, the variables <code>dofs_per_cell</code>
-// and <code>n_q_points</code> now are 9 each, where they were 4
-// before. Introducing such variables as abbreviations is a good strategy to
-// make code work with different elements without having to change too much
-// code. Secondly, the <code>fe_values</code> object of course needs to do
-// other things as well, since the shape functions are now quadratic, rather
-// than linear, in each coordinate variable. Again, however, this is something
-// that is completely transparent to user code and nothing that you have to
-// worry about.
-
-template <int dim>
-void Step6<dim>::assemble_system ()
-{
-  system_matrix = 0;
-  const QGauss<dim>  qformula(fe.degree+1);
-
-  FEValues<dim> fe_values (fe, qformula,
-                           update_values    |  update_gradients |
-                           update_quadrature_points  |  update_JxW_values);
-
-  FullMatrix<double>   cell_matrix;
-
-  Vector<double>       cell_rhs;
-
-
-  std::vector<unsigned int> local_dof_indices;
-
-
-  const Coefficient<dim> coefficient;
-  std::vector<double>    coefficient_values;
-
-  typename 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_matrix = 0;
-      cell_rhs.reinit (dofs_per_cell);
-      cell_rhs = 0;
-
-      fe_values.reinit (cell);
-
-      coefficient_values.resize(fe_values.n_quadrature_points);
-
-      coefficient.value_list (fe_values.get_quadrature_points(),
-                              coefficient_values);
-
-      for (unsigned int q_point=0; q_point<fe_values.n_quadrature_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) *
-                            1.0 *
-                            fe_values.JxW(q_point));
-          }
-
-      // Finally, transfer the contributions from @p cell_matrix and
-      // @p cell_rhs into the global objects.
-      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);
-    }
-  // Now we are done assembling the linear system.  The constrained nodes are
-  // still in the linear system (there is a one on the diagonal of the matrix
-  // and all other entries for this line are set to zero) but the computed
-  // values are invalid. We compute the correct values for these nodes at the
-  // end of the <code>solve</code> function.
-}
-
-
-template <int dim>
-void Step6<dim>::assemble_system_hp ()
-{
-  system_matrix = 0;
-
-  //  const QGauss<dim>  quadrature_formula(3);
-  hp::QCollection<dim> qformulas;
-  qformulas.push_back(QGauss<dim>(fe.degree+1));
-
-  hp::FEValues<dim> hp_fe_values (hpfe, qformulas,
-                                  update_values    |  update_gradients |
-                                  update_quadrature_points  |  update_JxW_values);
-
-  FullMatrix<double>   cell_matrix;
-
-  Vector<double>       cell_rhs;
-
-
-  std::vector<unsigned int> local_dof_indices;
-
-
-  const Coefficient<dim> coefficient;
-  std::vector<double>    coefficient_values;
-
-  typename hp::DoFHandler<dim>::active_cell_iterator
-  cell = hpdof_handler.begin_active(),
-  endc = hpdof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      hp_fe_values.reinit (cell);
-      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
-
-      const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
-      cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
-      cell_matrix = 0;
-      cell_rhs.reinit (dofs_per_cell);
-      cell_rhs = 0;
-
-      //  fe_values.reinit (cell);
-
-      coefficient_values.resize(fe_values.n_quadrature_points);
-
-      coefficient.value_list (fe_values.get_quadrature_points(),
-                              coefficient_values);
-
-      for (unsigned int q_point=0; q_point<fe_values.n_quadrature_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) *
-                            1.0 *
-                            fe_values.JxW(q_point));
-          }
-
-      // Finally, transfer the contributions from @p cell_matrix and
-      // @p cell_rhs into the global objects.
-      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);
-    }
-  // Now we are done assembling the linear system.  The constrained nodes are
-  // still in the linear system (there is a one on the diagonal of the matrix
-  // and all other entries for this line are set to zero) but the computed
-  // values are invalid. We compute the correct values for these nodes at the
-  // end of the <code>solve</code> function.
-}
-
-
-
-// @sect4{Step6::solve}
-
-// We continue with gradual improvements. The function that solves the linear
-// system again uses the SSOR preconditioner, and is again unchanged except
-// that we have to incorporate hanging node constraints. As mentioned above,
-// the degrees of freedom from the ConstraintMatrix corresponding to hanging
-// node constraints and boundary values have been removed from the linear
-// system by giving the rows and columns of the matrix a special
-// treatment. This way, the values for these degrees of freedom have wrong,
-// but well-defined values after solving the linear system. What we then have
-// to do is to use the constraints to assign to them the values that they
-// should have. This process, called <code>distributing</code> constraints,
-// computes the values of constrained nodes from the values of the
-// unconstrained ones, and requires only a single additional function call
-// that you find at the end of this function:
-
-template <int dim>
-void Step6<dim>::solve ()
-{
-  /*
-  SparseDirectUMFPACK u;
-  u.initialize(system_matrix);
-  u.vmult(solution, system_rhs);
-  */
-  /*
-
-
-
-  SolverControl           solver_control (10000, 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);
-}
-
-
-// @sect4{Step6::refine_grid}
-
-// Instead of global refinement, we now use a slightly more elaborate
-// scheme. We will use the <code>KellyErrorEstimator</code> class which
-// implements an error estimator for the Laplace equation; it can in principle
-// handle variable coefficients, but we will not use these advanced features,
-// but rather use its most simple form since we are not interested in
-// quantitative results but only in a quick way to generate locally refined
-// grids.
-//
-// Although the error estimator derived by Kelly et al. was originally
-// developed for the Laplace equation, we have found that it is also well
-// suited to quickly generate locally refined grids for a wide class of
-// problems. Basically, it looks at the jumps of the gradients of the solution
-// over the faces of cells (which is a measure for the second derivatives) and
-// scales it by the size of the cell. It is therefore a measure for the local
-// smoothness of the solution at the place of each cell and it is thus
-// understandable that it yields reasonable grids also for hyperbolic
-// transport problems or the wave equation as well, although these grids are
-// certainly suboptimal compared to approaches specially tailored to the
-// problem. This error estimator may therefore be understood as a quick way to
-// test an adaptive program.
-//
-// The way the estimator works is to take a <code>DoFHandler</code> object
-// describing the degrees of freedom and a vector of values for each degree of
-// freedom as input and compute a single indicator value for each active cell
-// of the triangulation (i.e. one value for each of the
-// <code>triangulation.n_active_cells()</code> cells). To do so, it needs two
-// additional pieces of information: a quadrature formula on the faces
-// (i.e. quadrature formula on <code>dim-1</code> dimensional objects. We use
-// a 3-point Gauss rule again, a pick that is consistent and appropriate with
-// the choice bi-quadratic finite element shape functions in this program.
-// (What constitutes a suitable quadrature rule here of course depends on
-// knowledge of the way the error estimator evaluates the solution field. As
-// said above, the jump of the gradient is integrated over each face, which
-// would be a quadratic function on each face for the quadratic elements in
-// use in this example. In fact, however, it is the square of the jump of the
-// gradient, as explained in the documentation of that class, and that is a
-// quartic function, for which a 3 point Gauss formula is sufficient since it
-// integrates polynomials up to order 5 exactly.)
-//
-// Secondly, the function wants a list of boundaries where we have imposed
-// Neumann value, and the corresponding Neumann values. This information is
-// represented by an object of type <code>FunctionMap::type</code> that is
-// essentially a map from boundary indicators to function objects describing
-// Neumann boundary values (in the present example program, we do not use
-// Neumann boundary values, so this map is empty, and in fact constructed
-// using the default constructor of the map in the place where the function
-// call expects the respective function argument).
-//
-// The output, as mentioned is a vector of values for all cells. While it may
-// make sense to compute the *value* of a degree of freedom very accurately,
-// it is usually not helpful to compute the *error indicator* corresponding to
-// a cell particularly accurately. We therefore typically use a vector of
-// floats instead of a vector of doubles to represent error indicators.
-template <int dim>
-void Step6<dim>::refine_grid ()
-{
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-  for (unsigned int i=0; i<triangulation.n_active_cells(); ++i)
-    estimated_error_per_cell (i)=i*1.0;
-
-  /*
-  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.0);
-
-  // After the previous function has exited, some cells are flagged for
-  // refinement, and some other for coarsening. The refinement or coarsening
-  // itself is not performed by now, however, since there are cases where
-  // further modifications of these flags is useful. Here, we don't want to do
-  // any such thing, so we can tell the triangulation to perform the actions
-  // for which the cells are flagged:
-  triangulation.execute_coarsening_and_refinement ();
-}
-
-
-// @sect4{Step6::output_results}
-
-// At the end of computations on each grid, and just before we continue the
-// next cycle with mesh refinement, we want to output the results from this
-// cycle.
-//
-// In the present program, we will not write the solution (except for in the
-// last step, see the next function), but only the meshes that we generated,
-// as a two-dimensional Encapsulated Postscript (EPS) file.
-//
-// We have already seen in step-1 how this can be achieved. The only thing we
-// have to change is the generation of the file name, since it should contain
-// the number of the present refinement cycle provided to this function as an
-// argument. The most general way is to use the std::stringstream class as
-// shown in step-5, but here's a little hack that makes it simpler if we know
-// that we have less than 10 iterations: assume that the %numbers `0' through
-// `9' are represented consecutively in the character set used on your machine
-// (this is in fact the case in all known character sets), then '0'+cycle
-// gives the character corresponding to the present cycle number. Of course,
-// this will only work if the number of cycles is actually less than 10, and
-// rather than waiting for the disaster to happen, we safeguard our little
-// hack with an explicit assertion at the beginning of the function. If this
-// assertion is triggered, i.e. when <code>cycle</code> is larger than or
-// equal to 10, an exception of type <code>ExcNotImplemented</code> is raised,
-// indicating that some functionality is not implemented for this case (the
-// functionality that is missing, of course, is the generation of file names
-// for that case):
-template <int dim>
-void Step6<dim>::output_results (const unsigned int cycle) const
-{
-  Assert (cycle < 10, ExcNotImplemented());
-
-  std::string filename = "grid-";
-  filename += ('0' + cycle);
-  filename += ".eps";
-
-  std::ofstream output (filename.c_str());
-
-  GridOut grid_out;
-  grid_out.write_eps (triangulation, output);
-}
-
-
-
-// @sect4{Step6::run}
-
-// The final function before <code>main()</code> is again the main driver of
-// the class, <code>run()</code>. It is similar to the one of step-5, except
-// that we generate a file in the program again instead of reading it from
-// disk, in that we adaptively instead of globally refine the mesh, and that
-// we output the solution on the final mesh in the present function.
-//
-// The first block in the main loop of the function deals with mesh
-// generation. If this is the first cycle of the program, instead of reading
-// the grid from a file on disk as in the previous example, we now again
-// create it using a library function. The domain is again a circle, which is
-// why we have to provide a suitable boundary object as well. We place the
-// center of the circle at the origin and have the radius be one (these are
-// the two hidden arguments to the function, which have default values).
-//
-// You will notice by looking at the coarse grid that it is of inferior
-// quality than the one which we read from the file in the previous example:
-// the cells are less equally formed. However, using the library function this
-// program works in any space dimension, which was not the case before.
-//
-// In case we find that this is not the first cycle, we want to refine the
-// grid. Unlike the global refinement employed in the last example program, we
-// now use the adaptive procedure described above.
-//
-// The rest of the loop looks as before:
-template <int dim>
-void Step6<dim>::run ()
-{
-  for (unsigned int cycle=0; cycle<2; ++cycle)
-    {
-      std::cout << "Cycle " << cycle << ':' << std::endl;
-
-      if (cycle == 0)
-        {
-          GridGenerator::hyper_ball (triangulation);
-
-          static const HyperBallBoundary<dim> boundary;
-          triangulation.set_boundary (0, boundary);
-
-          triangulation.refine_global (8);
-        }
-      else
-        refine_grid ();
-
-      std::cout << "   Number of active cells:       "
-                << triangulation.n_active_cells()
-                << std::endl;
-
-      std::cout << "setup" << std::endl;
-
-//    computing_timer.enter_section ("setup");
-      setup_system ();
-//    computing_timer.exit_section ("setup");
-
-      std::cout << "   Number of degrees of freedom: "
-                << dof_handler.n_dofs() << " " << hpdof_handler.n_dofs()
-                << std::endl;
-
-      std::cout << "warmup" << std::endl;
-      assemble_system ();// warm up
-      computing_timer.enter_section ("assembly");
-
-      std::cout << "assemble" << std::endl;
-
-      assemble_system ();
-      computing_timer.exit_section ("assembly");
-
-      std::cout << "assemble hp" << std::endl;
-
-      assemble_system_hp (); //warm up
-      computing_timer.enter_section ("assembly_hp");
-
-      assemble_system_hp ();
-      computing_timer.exit_section ("assembly_hp");
-
-      std::cout << "solve" << std::endl;
-
-      solve ();
-      //      output_results (cycle);
-      std::cout << "done" << std::endl;
-    }
-
-  // After we have finished computing the solution on the finest mesh, and
-  // writing all the grids to disk, we want to also write the actual solution
-  // on this final mesh to a file. As already done in one of the previous
-  // examples, we use the EPS format for output, and to obtain a reasonable
-  // view on the solution, we rescale the z-axis by a factor of four.
-  /*
-  DataOutBase::EpsFlags eps_flags;
-  eps_flags.z_scaling = 4;
-
-  DataOut<dim,hp::DoFHandler<dim> > data_out;
-  data_out.set_flags (eps_flags);
-
-  data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (solution, "solution");
-  data_out.build_patches ();
-
-  std::ofstream output ("final-solution.eps");
-  data_out.write_eps (output);*/
-}
-
-
-// @sect3{The <code>main</code> function}
-
-// The main function is unaltered in its functionality from the previous
-// example, but we have taken a step of additional caution. Sometimes,
-// something goes wrong (such as insufficient disk space upon writing an
-// output file, not enough memory when trying to allocate a vector or a
-// matrix, or if we can't read from or write to a file for whatever reason),
-// and in these cases the library will throw exceptions. Since these are
-// run-time problems, not programming errors that can be fixed once and for
-// all, this kind of exceptions is not switched off in optimized mode, in
-// contrast to the <code>Assert</code> macro which we have used to test
-// against programming errors. If uncaught, these exceptions propagate the
-// call tree up to the <code>main</code> function, and if they are not caught
-// there either, the program is aborted. In many cases, like if there is not
-// enough memory or disk space, we can't do anything but we can at least print
-// some text trying to explain the reason why the program failed. A way to do
-// so is shown in the following. It is certainly useful to write any larger
-// program in this way, and you can do so by more or less copying this
-// function except for the <code>try</code> block that actually encodes the
-// functionality particular to the present application.
-int main ()
-{
-
-  // The general idea behind the layout of this function is as follows: let's
-  // try to run the program as we did before...
-  try
-    {
-      deallog.depth_console (0);
-
-      Step6<2> laplace_problem_2d;
-      laplace_problem_2d.run ();
-    }
-  // ...and if this should fail, try to gather as much information as
-  // possible. Specifically, if the exception that was thrown is an object of
-  // a class that is derived from the C++ standard class
-  // <code>exception</code>, then we can use the <code>what</code> member
-  // function to get a string which describes the reason why the exception was
-  // thrown.
-  //
-  // The deal.II exception classes are all derived from the standard class,
-  // and in particular, the <code>exc.what()</code> function will return
-  // approximately the same string as would be generated if the exception was
-  // thrown using the <code>Assert</code> macro. You have seen the output of
-  // such an exception in the previous example, and you then know that it
-  // contains the file and line number of where the exception occured, and
-  // some other information. This is also what the following statements would
-  // print.
-  //
-  // Apart from this, there isn't much that we can do except exiting the
-  // program with an error code (this is what the <code>return 1;</code>
-  // does):
-  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;
-    }
-  // If the exception that was thrown somewhere was not an object of a class
-  // derived from the standard <code>exception</code> class, then we can't do
-  // anything at all. We then simply print an error message and exit.
-  catch (...)
-    {
-      std::cerr << std::endl << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-                << "Aborting!" << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      return 1;
-    }
-
-  // If we got to this point, there was no exception which propagated up to
-  // the main function (there may have been exceptions, but they were caught
-  // somewhere in the program or the library). Therefore, the program
-  // performed as was expected and we can return without error.
-  return 0;
-}
diff --git a/tests/benchmarks/test_poisson/CMakeLists.txt b/tests/benchmarks/test_poisson/CMakeLists.txt
deleted file mode 100644 (file)
index bc4ea38..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-##
-#  CMake script for the step-1 tutorial program:
-##
-
-# Set the name of the project and target:
-SET(TARGET "poisson")
-
-# Declare all source files the target consists of:
-SET(TARGET_SRC
-  ${TARGET}.cc
-  # You can specify additional files here!
-  )
-
-# Usually, you will not need to modify anything beyond this point...
-
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-FIND_PACKAGE(deal.II 8.0 QUIET
-  HINTS
-    ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR}
-  #
-  # If the deal.II library cannot be found (because it is not installed at a
-  # default location or your project resides at an uncommon place), you
-  # can specify additional hints for search paths here, e.g.
-  # "$ENV{HOME}/workspace/deal.II"
-  )
-
-IF (NOT ${deal.II_FOUND})
-   MESSAGE(FATAL_ERROR
-           "\n\n"
-          " *** Could not locate deal.II. *** "
-          "\n\n"
-           " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag 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/tests/benchmarks/test_poisson/poisson.cc b/tests/benchmarks/test_poisson/poisson.cc
deleted file mode 100644 (file)
index 959809b..0000000
+++ /dev/null
@@ -1,242 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2013 - 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.
-//
-// ---------------------------------------------------------------------
-
-
-
-const unsigned int element_degree = 2;
-const unsigned int dimension = 3;
-
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/timer.h>
-#include <deal.II/base/function.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/compressed_simple_sparsity_pattern.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/dofs/dof_handler.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>
-
-
-using namespace dealii;
-
-
-template <int dim>
-class HelmholtzProblem
-{
-public:
-  HelmholtzProblem (const FiniteElement<dim> &fe);
-  void run ();
-
-private:
-  void setup_system ();
-  void assemble_system ();
-  void solve ();
-
-  Triangulation<dim>                      triangulation;
-  const FiniteElement<dim>               &fe;
-  DoFHandler<dim>                         dof_handler;
-
-  ConstraintMatrix                        hanging_node_constraints;
-
-  SparsityPattern                         sparsity_pattern;
-  SparseMatrix<double>                    system_matrix;
-  Vector<double>                          tri_sol, tri_rhs;
-  TimerOutput                             timer;
-};
-
-
-
-template <int dim>
-HelmholtzProblem<dim>::HelmholtzProblem (const FiniteElement<dim> &fe) :
-  fe (fe),
-  dof_handler (triangulation),
-  timer(std::cout, TimerOutput::summary, TimerOutput::wall_times)
-{}
-
-
-
-template <int dim>
-void HelmholtzProblem<dim>::setup_system ()
-{
-  timer.enter_subsection("setup mesh and matrix");
-
-  GridGenerator::hyper_cube (triangulation, -1, 1);
-  triangulation.refine_global(6);
-  dof_handler.distribute_dofs (fe);
-  std::cout << "Number of active cells:          "
-            << triangulation.n_active_cells()
-            << std::endl
-            << "Number total degrees of freedom: "
-            << dof_handler.n_dofs() << std::endl;
-
-  hanging_node_constraints.clear ();
-  IndexSet locally_relevant (dof_handler.locally_owned_dofs().size());
-  DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant);
-  hanging_node_constraints.reinit (locally_relevant);
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                           hanging_node_constraints);
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                            0,
-                                            ZeroFunction<dim>(),
-                                            hanging_node_constraints);
-
-  {
-    CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
-                                         dof_handler.n_dofs(),
-                                         locally_relevant);
-    DoFTools::make_sparsity_pattern (dof_handler, csp,
-                                     hanging_node_constraints, false);
-    sparsity_pattern.copy_from (csp);
-  }
-  system_matrix.reinit(sparsity_pattern);
-  tri_sol.reinit (dof_handler.n_dofs());
-  tri_rhs.reinit (tri_sol);
-
-  timer.leave_subsection();
-}
-
-
-template <int dim>
-void HelmholtzProblem<dim>::assemble_system ()
-{
-  timer.enter_subsection("write into matrix");
-
-  QGauss<dim>   quadrature_formula(fe.degree+1);
-
-  const unsigned int n_q_points    = quadrature_formula.size();
-  const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
-  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);
-
-  FEValues<dim>  fe_values (fe, quadrature_formula,
-                            update_values   | update_gradients |
-                            update_JxW_values);
-
-  typename DoFHandler<dim>::active_cell_iterator
-  cell = dof_handler.begin_active(),
-  endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      fe_values.reinit (cell);
-      if (fe_values.get_cell_similarity() != CellSimilarity::translation)
-        cell_matrix = 0;
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          {
-            if (fe_values.get_cell_similarity() != CellSimilarity::translation)
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
-                                      fe_values.shape_grad(j,q_point)
-                                      +
-                                      fe_values.shape_value(i,q_point) *
-                                      fe_values.shape_value(j,q_point)) *
-                                     fe_values.JxW(q_point));
-          }
-
-      cell->get_dof_indices (local_dof_indices);
-      hanging_node_constraints.distribute_local_to_global (cell_matrix,
-                                                           local_dof_indices,
-                                                           system_matrix);
-    }
-
-  timer.leave_subsection();
-}
-
-
-
-template <int dim>
-void HelmholtzProblem<dim>::solve ()
-{
-  for (unsigned int i=0; i<tri_rhs.size(); i++)
-    if (hanging_node_constraints.is_constrained(i)==false)
-      {
-        const double value = (double)myrand()/RAND_MAX;
-        tri_rhs(i) = value;
-      }
-
-  timer.enter_subsection("40 matrix-vector products");
-
-  for (unsigned int i=0; i<40; ++i)
-    {
-      system_matrix.vmult (tri_sol, tri_rhs);
-    }
-
-  timer.leave_subsection();
-}
-
-
-template <int dim>
-void HelmholtzProblem<dim>::run ()
-{
-  setup_system();
-  assemble_system();
-  solve();
-}
-
-int main (int argc, char **argv)
-{
-
-  try
-    {
-      Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());
-      deallog.depth_console (0);
-
-      FE_Q<dimension> fe(element_degree);
-      HelmholtzProblem<dimension> problem(fe);
-      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;
-      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/tests/benchmarks/testlist.sh b/tests/benchmarks/testlist.sh
deleted file mode 100755 (executable)
index 96858ae..0000000
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash
-export TESTS="step-22 tablehandler test_assembly test_poisson test_hp"

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.