+++ /dev/null
-#!/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"
-
+++ /dev/null
-#!/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`"
+++ /dev/null
-#!/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
+++ /dev/null
-#!/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
+++ /dev/null
-
-gettimes: get_times.cc Makefile
- @g++ get_times.cc -o gettimes
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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;
-}
+++ /dev/null
-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
+++ /dev/null
-#!/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"
-
+++ /dev/null
-#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`"
-
+++ /dev/null
-#!/bin/bash
-
-
-while :
-do
- ./bench.sh
-
- ./doplots.sh
-
- ./sync.sh
-
- sleep 3600
-done
\ No newline at end of file
+++ /dev/null
-#!/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
-
+++ /dev/null
-##
-# 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()
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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 – 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 — 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;
-}
+++ /dev/null
-#!/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
-
+++ /dev/null
-##
-# 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()
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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");
-}
+++ /dev/null
-##
-# 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()
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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 – 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;
-}
+++ /dev/null
-##
-# 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()
+++ /dev/null
-/* ---------------------------------------------------------------------
- *
- * 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;
-}
+++ /dev/null
-##
-# 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()
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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;
-}
-
-
-
+++ /dev/null
-#!/bin/bash
-export TESTS="step-22 tablehandler test_assembly test_poisson test_hp"