From 7a9808e94b9dcbdbbec249e713516c6a3424d840 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Thu, 7 Jul 2016 12:19:28 +0100 Subject: [PATCH] remove tests/benchmarks/ they have been located at https://github.com/dealii/performance- benchmarks for a long time --- tests/benchmarks/baselineplot.sh | 36 - tests/benchmarks/bench.sh | 28 - tests/benchmarks/benchrev.sh | 40 - tests/benchmarks/doplots.sh | 20 - tests/benchmarks/gettimes/Makefile | 3 - tests/benchmarks/gettimes/get_times.cc | 166 -- tests/benchmarks/interactive.py | 155 -- tests/benchmarks/plot.sh | 28 - tests/benchmarks/qsub.regression | 88 -- tests/benchmarks/run.sh | 13 - tests/benchmarks/setup.sh | 28 - tests/benchmarks/step-22/CMakeLists.txt | 39 - tests/benchmarks/step-22/step-22.cc | 1386 ----------------- tests/benchmarks/sync.sh | 6 - tests/benchmarks/tablehandler/CMakeLists.txt | 39 - .../benchmarks/tablehandler/table_handler.cc | 63 - tests/benchmarks/test_assembly/CMakeLists.txt | 39 - tests/benchmarks/test_assembly/step-22.cc | 1339 ---------------- tests/benchmarks/test_hp/CMakeLists.txt | 39 - tests/benchmarks/test_hp/step.cc | 931 ----------- tests/benchmarks/test_poisson/CMakeLists.txt | 39 - tests/benchmarks/test_poisson/poisson.cc | 242 --- tests/benchmarks/testlist.sh | 2 - 23 files changed, 4769 deletions(-) delete mode 100755 tests/benchmarks/baselineplot.sh delete mode 100755 tests/benchmarks/bench.sh delete mode 100755 tests/benchmarks/benchrev.sh delete mode 100755 tests/benchmarks/doplots.sh delete mode 100644 tests/benchmarks/gettimes/Makefile delete mode 100644 tests/benchmarks/gettimes/get_times.cc delete mode 100755 tests/benchmarks/interactive.py delete mode 100755 tests/benchmarks/plot.sh delete mode 100644 tests/benchmarks/qsub.regression delete mode 100755 tests/benchmarks/run.sh delete mode 100755 tests/benchmarks/setup.sh delete mode 100644 tests/benchmarks/step-22/CMakeLists.txt delete mode 100644 tests/benchmarks/step-22/step-22.cc delete mode 100755 tests/benchmarks/sync.sh delete mode 100644 tests/benchmarks/tablehandler/CMakeLists.txt delete mode 100644 tests/benchmarks/tablehandler/table_handler.cc delete mode 100644 tests/benchmarks/test_assembly/CMakeLists.txt delete mode 100644 tests/benchmarks/test_assembly/step-22.cc delete mode 100644 tests/benchmarks/test_hp/CMakeLists.txt delete mode 100644 tests/benchmarks/test_hp/step.cc delete mode 100644 tests/benchmarks/test_poisson/CMakeLists.txt delete mode 100644 tests/benchmarks/test_poisson/poisson.cc delete mode 100755 tests/benchmarks/testlist.sh diff --git a/tests/benchmarks/baselineplot.sh b/tests/benchmarks/baselineplot.sh deleted file mode 100755 index 8ae6bd3191..0000000000 --- a/tests/benchmarks/baselineplot.sh +++ /dev/null @@ -1,36 +0,0 @@ -#!/bin/bash - - -source testlist.sh - - -cat </dev/null - echo -n "" > temp.txt - for a in {1..5}; do - echo "*" >> temp.txt - make run | grep "|" >> temp.txt - done - ./../gettimes/gettimes > names.test - if [[ -s names.test ]] ; then - words=`wc -w names.test | cut -f1 -d' '` - if [ "$words" -gt "0" ] ; then - cp names.test ../names.$test - fi ; - rm -rf names.test - fi ; - ./../gettimes/gettimes $PREVREVISION >>../datatable.$test - cd .. - - done diff --git a/tests/benchmarks/doplots.sh b/tests/benchmarks/doplots.sh deleted file mode 100755 index 7358c0456a..0000000000 --- a/tests/benchmarks/doplots.sh +++ /dev/null @@ -1,20 +0,0 @@ -#!/bin/bash - -source testlist.sh - -echo "generating images..." - -for test in $TESTS ; do - LASTREV=`tail -n 1 datatable.$test | cut -f 1 -d ' '` - ./plot.sh $test $LASTREV >script - gnuplot script - rm -rf script - convert -density 150 $test.eps $test.png -done - - -./baselineplot.sh > script -gnuplot script -rm -rf script -convert -density 150 baseline.eps baseline.png -python interactive.py >index.html \ No newline at end of file diff --git a/tests/benchmarks/gettimes/Makefile b/tests/benchmarks/gettimes/Makefile deleted file mode 100644 index 3efb00079e..0000000000 --- a/tests/benchmarks/gettimes/Makefile +++ /dev/null @@ -1,3 +0,0 @@ - -gettimes: get_times.cc Makefile - @g++ get_times.cc -o gettimes diff --git a/tests/benchmarks/gettimes/get_times.cc b/tests/benchmarks/gettimes/get_times.cc deleted file mode 100644 index 35d4ab89e6..0000000000 --- a/tests/benchmarks/gettimes/get_times.cc +++ /dev/null @@ -1,166 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2013 - 2014 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -#include -#include -#include -#include -#include -#include -#include - -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 names; - vector times; - - //If no revision number, retrieve column names - if (argc <= 1) - { - string curr_line; - - //Ignore first n lines + delimeter- they don't hold any useful data - for (int i = 0; i < IGNORE_LINES + 1; i++) - { - getline(input, curr_line); - } - - while (!input.eof()) - { - - getline(input, curr_line); - if (curr_line == "") - continue; - if (curr_line == DELIM) - break; - - ////cout << "curr line: " << curr_line << endl; - //Looking for the string after the '|' and ' ' - int first_char = 0; - if (curr_line[0] == '|') - first_char++; - if (curr_line[1] == ' ') - first_char++; - - - ////cout << "first char at pos: " << first_char << endl; - - //Find end of string - int num_chars = 0; - while (curr_line[first_char + num_chars] != '|') - num_chars++; - - num_chars--; - - while (curr_line[first_char + num_chars] == ' ') - num_chars--; - - names.push_back(curr_line.substr(first_char, num_chars+1)); - } - } - - else //Else, extract execution time from each line - { - - int time_index = 0; - - while (!input.eof()) - { - - string curr_line = ""; - - //Read in line - getline(input,curr_line); - - //Check for delimeter - if (curr_line == DELIM) - { - ////cout << "Delimeter detected" << endl; - time_index = 0; - //Ignore first n lines- they don't hold any useful data - for (int i = 0; i < IGNORE_LINES; i++) - { - string dummy; - getline(input, dummy); - ////cout << "Skipping: " << dummy << endl; - } - continue; - } - - ////cout << "Reading: " << curr_line << endl; - - if (curr_line == "" && DELIM != "") - continue; - - - //Looking for a number that ends with 's' - int last_s = curr_line.rfind('s'); - //In case 's' is not used in the future - assert(isdigit(curr_line[last_s - 1])); // Test: s preceded by number - - //Find time string and convert to double - int num_start = last_s - 1; - while (curr_line[num_start] != ' ') - num_start--; - string timestr = curr_line.substr(num_start+1, last_s - num_start); - double time = (double)atof(timestr.substr(0,timestr.size()-1).c_str()); - - assert(times.size() >= time_index); - // first addition of times to vector; each loop, times.size() should be one less than time_index - if (times.size() == time_index) - times.push_back(time); - else - times[time_index] = (time < times[time_index]) ? time : times[time_index]; //else, determines minimum time and stores it - time_index++; - - } - } - - - //Output individual names - if (argc <= 1) - { - for (int i = 0; i < names.size(); i++) - { - cout << names[i]; - if (i < names.size() - 1) - cout << endl; - } - } - - //Output individual times - if (argc > 1) - { - cout << argv[1]; - for (int i = 0; i < times.size(); i++) - cout << " " << times[i]; - } - - cout << endl; - - - return 0; -} diff --git a/tests/benchmarks/interactive.py b/tests/benchmarks/interactive.py deleted file mode 100755 index f687f4b079..0000000000 --- a/tests/benchmarks/interactive.py +++ /dev/null @@ -1,155 +0,0 @@ -import textwrap -import os - -begin = \ -""" - - - - - deal.II regression timings - - - - - - - - - -deal.II performance benchmarks, see -http://www.dealii.org/testsuite.html

- -
- - - - - - - - - - - - - - -""" - - - -list = os.listdir(".") - -print begin - -first = 1 -for fname in list: - if (fname.startswith("names.")): - testname = fname[6:] - names = open(fname).readlines() - data = open("datatable."+testname).readlines() - idx = 0 - for name in names: - if first == 1: - first = 0 - else: - print "," - idx = idx+1 - print "{ name: '%s - %s', data: [" % (testname,name[:-1]) - i=0 -#((\$$col-$baseline)/$baseline*100.0) - baseline = -1 - for l in data: - if (len(l.strip())<1): - continue - - if (baseline>-1 and baseline<0.5): - continue; - - def isfloat(x): - try: - float(x) - return True - except ValueError: - return False - lnumbers = [float(x) for x in l.split() if isfloat(x)] - if len(lnumbers)<=1: - continue; - - if baseline<0 and len(lnumbers)>idx: - baseline=lnumbers[idx] - - if lnumbers[0]<27000: - continue; - - if (i==1): - print(","), - i=1; - - - if len(lnumbers)>idx: - print "[%d,%f]" % (lnumbers[0], (lnumbers[idx]-baseline)/baseline*100.0) - print "]}\n" - - - - -print end diff --git a/tests/benchmarks/plot.sh b/tests/benchmarks/plot.sh deleted file mode 100755 index ec9b45278c..0000000000 --- a/tests/benchmarks/plot.sh +++ /dev/null @@ -1,28 +0,0 @@ -#!/bin/bash - -#launch with the name of test to generate the .eps for - -cat < -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -// Then we need to include the header file -// for the sparse direct solver UMFPACK: -#include -#include - -// This includes the library for the -// incomplete LU factorization that will -// be used as a preconditioner in 3D: -#include - -// This is C++: -#include -#include - -// 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 - struct InnerPreconditioner; - - // In 2D, we are going to use a sparse direct - // solver as preconditioner: - template <> - struct InnerPreconditioner<2> - { - typedef SparseILU type; -// typedef SparseDirectUMFPACK type; - }; - - // And the ILU preconditioning in 3D, called - // by SparseILU: - template <> - struct InnerPreconditioner<3> - { - typedef SparseILU type; - }; - - - // @sect3{The StokesProblem 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 - // hanging_node_constraints - // into constraints. - template - 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 triangulation; - FESystem fe; - DoFHandler dof_handler; - - ConstraintMatrix constraints; - - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; - - BlockVector solution; - BlockVector 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::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 - class BoundaryValues : public Function - { - public: - BoundaryValues () : Function(dim+1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; - }; - - - template - double - BoundaryValues::value (const Point &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 - void - BoundaryValues::vector_value (const Point &p, - Vector &values) const - { - for (unsigned int c=0; cn_components; ++c) - values(c) = BoundaryValues::value (p, c); - } - - - - // We implement similar functions for - // the right hand side which for the - // current example is simply zero: - template - class RightHandSide : public Function - { - public: - RightHandSide () : Function(dim+1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; - - }; - - - template - double - RightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const - { - return 0; - } - - - template - void - RightHandSide::vector_value (const Point &p, - Vector &values) const - { - for (unsigned int c=0; cn_components; ++c) - values(c) = RightHandSide::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 InverseMatrix class template} - - // The InverseMatrix - // 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 - // InverseMatrix object - // is created. The member function - // vmult is, as in - // step-20, a multiplication with a - // vector, obtained by solving a - // linear system: - template - class InverseMatrix : public Subscriptor - { - public: - InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer matrix; - const SmartPointer preconditioner; - }; - - - template - InverseMatrix::InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner) - : - matrix (&m), - preconditioner (&preconditioner) - {} - - - // This is the implementation of the - // vmult 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 - void InverseMatrix::vmult (Vector &dst, - const Vector &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 SchurComplement 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 - // Preconditioner in order to - // access that when specifying the respective - // type of the inverse matrix class. As a - // consequence of the definition above, the - // declaration InverseMatrix now - // contains the second template parameter - // for a preconditioner class as above, which - // affects the SmartPointer - // object m_inverse as well. - template - class SchurComplement : public Subscriptor - { - public: - SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix, Preconditioner> &A_inverse); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer > system_matrix; - const SmartPointer, Preconditioner> > A_inverse; - - mutable Vector tmp1, tmp2; - }; - - - - template - SchurComplement:: - SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix,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 - void SchurComplement::vmult (Vector &dst, - const Vector &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 - // degree+1 for the - // vector-valued velocity components - // and of order degree - // 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 - // Triangulation::MeshSmoothing - // for details. - template - StokesProblem::StokesProblem (const unsigned int degree) - : - degree (degree), - triangulation (Triangulation::maximum_smoothing), - fe (FE_Q(degree+1), dim, - FE_Q(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 - // DoFRenumbering::component_wise - // 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 - // section on improved ILU - // 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 - // dim 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 - // block_component 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 - void StokesProblem::setup_dofs () - { - A_preconditioner.reset (); - system_matrix.clear (); - - dof_handler.distribute_dofs (fe); - DoFRenumbering::Cuthill_McKee (dof_handler); - - std::vector 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 component_mask 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 - // H1 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 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(), - 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 - // DoFTools::count_dofs_per_block - // does the same as - // DoFTools::count_dofs_per_component, - // but now grouped as velocity and pressure - // block via block_component. - std::vector 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 csp - // will be released once the information - // has been copied to - // sparsity_pattern. - { - 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 - void StokesProblem::assemble_system () - { - system_matrix=0; - system_rhs=0; - - QGauss quadrature_formula(degree+2); - - FEValues 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 local_matrix (dofs_per_cell, dofs_per_cell); - Vector local_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const RightHandSide right_hand_side; - std::vector > rhs_values (n_q_points, - Vector(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 dofs_per_cell - // different values per quadrature - // point, but we extract it - // dofs_per_cell*dofs_per_cell - // times from the FEValues object - - // for both the loop over - // i and the inner - // loop over j. 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 dofs_per_cell, - // but only up to i, - // the index of the outer loop. - std::vector > symgrad_phi_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); - - typename DoFHandler::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 phi_p[i] * - // phi_p[j] , 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 phi_p[i] * - // phi_p[j] 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; iget_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., - // block(0,0) in the - // system matrix. As mentioned - // above, this depends on the - // spatial dimension. Since the two - // classes described by the - // InnerPreconditioner::type - // 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::type>(new typename InnerPreconditioner::type()); - A_preconditioner->initialize (system_matrix.block(0,0), - typename InnerPreconditioner::type::AdditionalData()); - - } - - - - // @sect4{StokesProblem::solve} - - // After the discussion in the introduction - // and the definition of the respective - // classes above, the implementation of the - // solve function is rather - // straigt-forward and done in a similar way - // as in step-20. To start with, we need an - // object of the InverseMatrix - // 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 - // InnerPreconditioner::type. - template - void StokesProblem::solve () - { - const InverseMatrix, - typename InnerPreconditioner::type> - A_inverse (system_matrix.block(0,0), *A_preconditioner); - Vector 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 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::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 - // InverseMatrix, which - // implements the vmult - // 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 preconditioner; - preconditioner.initialize (system_matrix.block(1,1), - SparseILU::AdditionalData()); - - InverseMatrix,SparseILU > - 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: velocity to the - // dim components of velocity - // and pressure 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 - // DataComponentInterpretation - // namespace: as with the filename, we - // create a vector in which the first - // dim components refer to the - // velocities and are given the tag - // DataComponentInterpretation::component_is_part_of_vector; - // we finally push one tag - // DataComponentInterpretation::component_is_scalar - // to describe the grouping of the pressure - // variable. - - // The rest of the function is then - // the same as in step-20. - template - void - StokesProblem::output_results (const unsigned int refinement_cycle) const - { - std::vector solution_names (dim, "velocity"); - solution_names.push_back ("pressure"); - - std::vector - data_component_interpretation - (dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - - DataOut data_out; - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, solution_names, - DataOut::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 StokesProblem 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 - void - StokesProblem::refine_mesh () - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - std::vector component_mask (dim+1, false); - component_mask[dim] = true; - KellyErrorEstimator::estimate (dof_handler, - QGauss(degree+1), - typename FunctionMap::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 - void StokesProblem::run () - { - { - std::vector subdivisions (dim, 1); - subdivisions[0] = 4; - - const Point bottom_left = (dim == 2 ? - Point(-2,-1) : - Point(-2,0,-1)); - const Point top_right = (dim == 2 ? - Point(2,0) : - Point(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::active_cell_iterator - cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell) - for (unsigned int f=0; f::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 main function} - -// The main function is the same as in -// step-20. We pass the element degree as a -// parameter and choose the space dimension -// at the well-known template slot. -int main () -{ - try - { - using namespace dealii; - using namespace Step22; - - deallog.depth_console (0); - - StokesProblem<2> flow_problem(1); - flow_problem.run (); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} diff --git a/tests/benchmarks/sync.sh b/tests/benchmarks/sync.sh deleted file mode 100755 index 1bc58c4d23..0000000000 --- a/tests/benchmarks/sync.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/bash -echo "copying files to public_html/" -cp *png ~/public_html/bench/ -cp index.html ~/public_html/bench/ -chmod a+r ~/public_html/bench/*png ~/public_html/bench/*html - diff --git a/tests/benchmarks/tablehandler/CMakeLists.txt b/tests/benchmarks/tablehandler/CMakeLists.txt deleted file mode 100644 index 41f57c92a1..0000000000 --- a/tests/benchmarks/tablehandler/CMakeLists.txt +++ /dev/null @@ -1,39 +0,0 @@ -## -# CMake script for the step-1 tutorial program: -## - -# Set the name of the project and target: -SET(TARGET "table_handler") - -# Declare all source files the target consists of: -SET(TARGET_SRC - ${TARGET}.cc - # You can specify additional files here! - ) - -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) - -FIND_PACKAGE(deal.II 8.0 QUIET - HINTS - ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR} - # - # If the deal.II library cannot be found (because it is not installed at a - # default location or your project resides at an uncommon place), you - # can specify additional hints for search paths here, e.g. - # "$ENV{HOME}/workspace/deal.II" - ) - -IF (NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR - "\n\n" - " *** Could not locate deal.II. *** " - "\n\n" - " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag to cmake \n" - " *** or set an environment variable \"DEAL_II_DIR\" that contains this path.") -ENDIF () - -DEAL_II_INITIALIZE_CACHED_VARIABLES() -PROJECT(${TARGET}) -DEAL_II_INVOKE_AUTOPILOT() diff --git a/tests/benchmarks/tablehandler/table_handler.cc b/tests/benchmarks/tablehandler/table_handler.cc deleted file mode 100644 index b5aaa660e2..0000000000 --- a/tests/benchmarks/tablehandler/table_handler.cc +++ /dev/null @@ -1,63 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2010 - 2014 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -// we test the performance of writing a large table - -#include -#include -#include -#include - -using namespace dealii; - - -#include -#include -#include -#include -#include - - -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 -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -// Then we need to include the header file -// for the sparse direct solver UMFPACK: -#include -#include - -// This includes the library for the -// incomplete LU factorization that will -// be used as a preconditioner in 3D: -#include - -// This is C++: -#include -#include - -// 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 - struct InnerPreconditioner; - - // In 2D, we are going to use a sparse direct - // solver as preconditioner: - template <> - struct InnerPreconditioner<2> - { - typedef SparseILU type; -// typedef SparseDirectUMFPACK type; - }; - - // And the ILU preconditioning in 3D, called - // by SparseILU: - template <> - struct InnerPreconditioner<3> - { - typedef SparseILU type; - }; - - - // @sect3{The StokesProblem 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 - // hanging_node_constraints - // into constraints. - template - 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 triangulation; - FESystem fe; - DoFHandler dof_handler; - - ConstraintMatrix constraints; - - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; - - BlockVector solution; - BlockVector 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::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 - class BoundaryValues : public Function - { - public: - BoundaryValues () : Function(dim+1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; - }; - - - template - double - BoundaryValues::value (const Point &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 - void - BoundaryValues::vector_value (const Point &p, - Vector &values) const - { - for (unsigned int c=0; cn_components; ++c) - values(c) = BoundaryValues::value (p, c); - } - - - - // We implement similar functions for - // the right hand side which for the - // current example is simply zero: - template - class RightHandSide : public Function - { - public: - RightHandSide () : Function(dim+1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; - - }; - - - template - double - RightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const - { - return 0; - } - - - template - void - RightHandSide::vector_value (const Point &p, - Vector &values) const - { - for (unsigned int c=0; cn_components; ++c) - values(c) = RightHandSide::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 InverseMatrix class template} - - // The InverseMatrix - // 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 - // InverseMatrix object - // is created. The member function - // vmult is, as in - // step-20, a multiplication with a - // vector, obtained by solving a - // linear system: - template - class InverseMatrix : public Subscriptor - { - public: - InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer matrix; - const SmartPointer preconditioner; - }; - - - template - InverseMatrix::InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner) - : - matrix (&m), - preconditioner (&preconditioner) - {} - - - // This is the implementation of the - // vmult 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 - void InverseMatrix::vmult (Vector &dst, - const Vector &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 SchurComplement 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 - // Preconditioner in order to - // access that when specifying the respective - // type of the inverse matrix class. As a - // consequence of the definition above, the - // declaration InverseMatrix now - // contains the second template parameter - // for a preconditioner class as above, which - // affects the SmartPointer - // object m_inverse as well. - template - class SchurComplement : public Subscriptor - { - public: - SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix, Preconditioner> &A_inverse); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer > system_matrix; - const SmartPointer, Preconditioner> > A_inverse; - - mutable Vector tmp1, tmp2; - }; - - - - template - SchurComplement:: - SchurComplement (const BlockSparseMatrix &system_matrix, - const InverseMatrix,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 - void SchurComplement::vmult (Vector &dst, - const Vector &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 - // degree+1 for the - // vector-valued velocity components - // and of order degree - // 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 - // Triangulation::MeshSmoothing - // for details. - template - StokesProblem::StokesProblem (const unsigned int degree) - : - degree (degree), - triangulation (Triangulation::maximum_smoothing), - fe (FE_Q(degree+1), dim, - FE_Q(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 - // DoFRenumbering::component_wise - // 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 - // section on improved ILU - // 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 - // dim 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 - // block_component 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 - void StokesProblem::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 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 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(), - 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 - // DoFTools::count_dofs_per_block - // does the same as - // DoFTools::count_dofs_per_component, - // but now grouped as velocity and pressure - // block via block_component. - std::vector 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 dsp - // will be released once the information has - // been copied to - // sparsity_pattern. - { - 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 - void StokesProblem::assemble_system () - { - system_matrix=0; - system_rhs=0; - - QGauss quadrature_formula(degree+2); - - FEValues 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 local_matrix (dofs_per_cell, dofs_per_cell); - Vector local_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const RightHandSide right_hand_side; - std::vector > rhs_values (n_q_points, - Vector(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 dofs_per_cell - // different values per quadrature - // point, but we extract it - // dofs_per_cell*dofs_per_cell - // times from the FEValues object - - // for both the loop over - // i and the inner - // loop over j. 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 dofs_per_cell, - // but only up to i, - // the index of the outer loop. - std::vector > symgrad_phi_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); - - typename DoFHandler::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 phi_p[i] * - // phi_p[j] , 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 phi_p[i] * - // phi_p[j] 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; iget_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., - // block(0,0) in the - // system matrix. As mentioned - // above, this depends on the - // spatial dimension. Since the two - // classes described by the - // InnerPreconditioner::type - // 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::type>(new typename InnerPreconditioner::type()); - A_preconditioner->initialize (system_matrix.block(0,0), - typename InnerPreconditioner::type::AdditionalData()); - - } - - - - // @sect4{StokesProblem::solve} - - // After the discussion in the introduction - // and the definition of the respective - // classes above, the implementation of the - // solve function is rather - // straigt-forward and done in a similar way - // as in step-20. To start with, we need an - // object of the InverseMatrix - // 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 - // InnerPreconditioner::type. - template - void StokesProblem::solve () - { - const InverseMatrix, - typename InnerPreconditioner::type> - A_inverse (system_matrix.block(0,0), *A_preconditioner); - Vector 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 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::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 - // InverseMatrix, which - // implements the vmult - // 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 preconditioner; - preconditioner.initialize (system_matrix.block(1,1), - SparseILU::AdditionalData()); - - InverseMatrix,SparseILU > - 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: velocity to the - // dim components of velocity - // and pressure 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 - // DataComponentInterpretation - // namespace: as with the filename, we - // create a vector in which the first - // dim components refer to the - // velocities and are given the tag - // DataComponentInterpretation::component_is_part_of_vector; - // we finally push one tag - // DataComponentInterpretation::component_is_scalar - // to describe the grouping of the pressure - // variable. - - // The rest of the function is then - // the same as in step-20. - template - void - StokesProblem::output_results (const unsigned int refinement_cycle) const - { - std::vector solution_names (dim, "velocity"); - solution_names.push_back ("pressure"); - - std::vector - data_component_interpretation - (dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - - DataOut data_out; - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, solution_names, - DataOut::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 StokesProblem 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 - void - StokesProblem::refine_mesh () - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - srand (0); - for (unsigned int i=0; i component_mask (dim+1, false); - component_mask[dim] = true; - KellyErrorEstimator::estimate (dof_handler, - QGauss(degree+1), - typename FunctionMap::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 - void StokesProblem::run () - { - { - std::vector subdivisions (dim, 1); - subdivisions[0] = 4; - - const Point bottom_left = (dim == 2 ? - Point(-2,-1) : - Point(-2,0,-1)); - const Point top_right = (dim == 2 ? - Point(2,0) : - Point(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::active_cell_iterator - cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell) - for (unsigned int f=0; f::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 main function} - -// The main function is the same as in -// step-20. We pass the element degree as a -// parameter and choose the space dimension -// at the well-known template slot. -int main () -{ - try - { - using namespace dealii; - using namespace Step22; - - deallog.depth_console (0); - - StokesProblem<2> flow_problem(1); - flow_problem.run (); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} diff --git a/tests/benchmarks/test_hp/CMakeLists.txt b/tests/benchmarks/test_hp/CMakeLists.txt deleted file mode 100644 index 0d9a9b43a6..0000000000 --- a/tests/benchmarks/test_hp/CMakeLists.txt +++ /dev/null @@ -1,39 +0,0 @@ -## -# CMake script for the step-1 tutorial program: -## - -# Set the name of the project and target: -SET(TARGET "step") - -# Declare all source files the target consists of: -SET(TARGET_SRC - ${TARGET}.cc - # You can specify additional files here! - ) - -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) - -FIND_PACKAGE(deal.II 8.0 QUIET - HINTS - ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR} - # - # If the deal.II library cannot be found (because it is not installed at a - # default location or your project resides at an uncommon place), you - # can specify additional hints for search paths here, e.g. - # "$ENV{HOME}/workspace/deal.II" - ) - -IF (NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR - "\n\n" - " *** Could not locate deal.II. *** " - "\n\n" - " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag to cmake \n" - " *** or set an environment variable \"DEAL_II_DIR\" that contains this path.") -ENDIF () - -DEAL_II_INITIALIZE_CACHED_VARIABLES() -PROJECT(${TARGET}) -DEAL_II_INVOKE_AUTOPILOT() diff --git a/tests/benchmarks/test_hp/step.cc b/tests/benchmarks/test_hp/step.cc deleted file mode 100644 index 0e47d46ccb..0000000000 --- a/tests/benchmarks/test_hp/step.cc +++ /dev/null @@ -1,931 +0,0 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2000 - 2015 by the deal.II authors - * - * This file is part of the deal.II library. - * - * The deal.II library is free software; you can use it, redistribute - * it, and/or modify it under the terms of the GNU Lesser General - * Public License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE at - * the top level of the deal.II distribution. - * - * --------------------------------------------------------------------- - - * Author: Wolfgang Bangerth, University of Heidelberg, 2000 - */ - -// @sect3{Include files} - -// The first few files have already been covered in previous examples and will -// thus not be further commented on. -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include -#include - -// From the following include file we will import the declaration of -// H1-conforming finite element shape functions. This family of finite -// elements is called FE_Q, 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 -// 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 -// grid_in.h: -#include - - -// When using locally refined grids, we will get so-called hanging -// nodes. 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 - -// 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 - -// 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 - -// Finally, this is as in previous programs: -using namespace dealii; - - -// @sect3{The Step6 class template} - -// The main class is again almost unchanged. Two additions, however, are made: -// we have added the refine_grid 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 -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 triangulation; - - DoFHandler dof_handler; - hp::DoFHandler hpdof_handler; - FE_Q fe; - hp::FECollection 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 system_matrix; - - Vector solution; - Vector system_rhs; -}; - - - -// @sect3{Nonconstant coefficients} - -// The implementation of nonconstant coefficients is copied verbatim from -// step-5: - -template -class Coefficient : public Function -{ -public: - Coefficient () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component = 0) const; -}; - - - -template -double Coefficient::value (const Point &p, - const unsigned int) const -{ - if (p.square() < 0.5*0.5) - return 20; - else - return 1; -} - - - -template -void Coefficient::value_list (const std::vector > &points, - std::vector &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; iStep6 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 1 in all previous examples) by -// the desired polynomial degree (here 2): -template -Step6::Step6 () - : - dof_handler (triangulation), - hpdof_handler (triangulation), - fe (3), - computing_timer (std::cout, - TimerOutput::summary, - TimerOutput::wall_times) - -{ - hpfe.push_back (FE_Q(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 dof_handler object was -// defined before and not after the fe 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 dof_handler.distribute_dofs(), -// the dof_handler 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 dof_handler object is destroyed, it would be unwise if we -// would allow the finite element object to be deleted before the -// dof_handler 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 -// Subscriptor/SmartPointer 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 Step6 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 fe object is -// destructed. What needs to be done is to tell the dof_handler -// object to release its lock to the finite element. Of course, the -// dof_handler 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 DoFHandler class has -// a function clear 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 -Step6::~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 -// dof_Handler.max_couplings_between_dofs() call will take care -// of this, however: -template -void Step6::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 constraints 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 constraints 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(), - constraints); - - - // The next step is closing 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 - // close() 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 constraints by setting the argument - // keep_constrained_dofs 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 true - // 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 QGauss 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 -// ConstraintMatrix::distribute_local_to_global 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 dofs_per_cell -// and n_q_points 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 fe_values 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 -void Step6::assemble_system () -{ - system_matrix = 0; - const QGauss qformula(fe.degree+1); - - FEValues fe_values (fe, qformula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - FullMatrix cell_matrix; - - Vector cell_rhs; - - - std::vector local_dof_indices; - - - const Coefficient coefficient; - std::vector coefficient_values; - - typename DoFHandler::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_pointget_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 solve function. -} - - -template -void Step6::assemble_system_hp () -{ - system_matrix = 0; - - // const QGauss quadrature_formula(3); - hp::QCollection qformulas; - qformulas.push_back(QGauss(fe.degree+1)); - - hp::FEValues hp_fe_values (hpfe, qformulas, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - FullMatrix cell_matrix; - - Vector cell_rhs; - - - std::vector local_dof_indices; - - - const Coefficient coefficient; - std::vector coefficient_values; - - typename hp::DoFHandler::active_cell_iterator - cell = hpdof_handler.begin_active(), - endc = hpdof_handler.end(); - for (; cell!=endc; ++cell) - { - hp_fe_values.reinit (cell); - const FEValues &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_pointget_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 solve 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 distributing 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 -void Step6::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 KellyErrorEstimator 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 DoFHandler 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 -// triangulation.n_active_cells() cells). To do so, it needs two -// additional pieces of information: a quadrature formula on the faces -// (i.e. quadrature formula on dim-1 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 FunctionMap::type 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 -void Step6::refine_grid () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - for (unsigned int i=0; i::estimate (dof_handler, - QGauss(3), - typename FunctionMap::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 cycle is larger than or -// equal to 10, an exception of type ExcNotImplemented 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 -void Step6::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 main() is again the main driver of -// the class, run(). 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 -void Step6::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 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 > 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 main 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 Assert macro which we have used to test -// against programming errors. If uncaught, these exceptions propagate the -// call tree up to the main 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 try 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 - // exception, then we can use the what 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 exc.what() function will return - // approximately the same string as would be generated if the exception was - // thrown using the Assert 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 return 1; - // 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 exception class, then we can't do - // anything at all. We then simply print an error message and exit. - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - // If we got to this point, there was no exception which propagated up to - // the main function (there may have been exceptions, but they were caught - // somewhere in the program or the library). Therefore, the program - // performed as was expected and we can return without error. - return 0; -} diff --git a/tests/benchmarks/test_poisson/CMakeLists.txt b/tests/benchmarks/test_poisson/CMakeLists.txt deleted file mode 100644 index bc4ea382fa..0000000000 --- a/tests/benchmarks/test_poisson/CMakeLists.txt +++ /dev/null @@ -1,39 +0,0 @@ -## -# CMake script for the step-1 tutorial program: -## - -# Set the name of the project and target: -SET(TARGET "poisson") - -# Declare all source files the target consists of: -SET(TARGET_SRC - ${TARGET}.cc - # You can specify additional files here! - ) - -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) - -FIND_PACKAGE(deal.II 8.0 QUIET - HINTS - ${deal.II_DIR}/ ${DEAL_II_DIR}/ ../../installed/ ../ ../../ ../../../ ../../../../../ $ENV{DEAL_II_DIR} - # - # If the deal.II library cannot be found (because it is not installed at a - # default location or your project resides at an uncommon place), you - # can specify additional hints for search paths here, e.g. - # "$ENV{HOME}/workspace/deal.II" - ) - -IF (NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR - "\n\n" - " *** Could not locate deal.II. *** " - "\n\n" - " *** You may want to either pass the -DDEAL_II_DIR=/path/to/deal.II flag to cmake \n" - " *** or set an environment variable \"DEAL_II_DIR\" that contains this path.") -ENDIF () - -DEAL_II_INITIALIZE_CACHED_VARIABLES() -PROJECT(${TARGET}) -DEAL_II_INVOKE_AUTOPILOT() diff --git a/tests/benchmarks/test_poisson/poisson.cc b/tests/benchmarks/test_poisson/poisson.cc deleted file mode 100644 index 959809b327..0000000000 --- a/tests/benchmarks/test_poisson/poisson.cc +++ /dev/null @@ -1,242 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2013 - 2015 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - - -const unsigned int element_degree = 2; -const unsigned int dimension = 3; - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -using namespace dealii; - - -template -class HelmholtzProblem -{ -public: - HelmholtzProblem (const FiniteElement &fe); - void run (); - -private: - void setup_system (); - void assemble_system (); - void solve (); - - Triangulation triangulation; - const FiniteElement &fe; - DoFHandler dof_handler; - - ConstraintMatrix hanging_node_constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - Vector tri_sol, tri_rhs; - TimerOutput timer; -}; - - - -template -HelmholtzProblem::HelmholtzProblem (const FiniteElement &fe) : - fe (fe), - dof_handler (triangulation), - timer(std::cout, TimerOutput::summary, TimerOutput::wall_times) -{} - - - -template -void HelmholtzProblem::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(), - 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 -void HelmholtzProblem::assemble_system () -{ - timer.enter_subsection("write into matrix"); - - QGauss 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 cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_JxW_values); - - typename DoFHandler::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_pointget_dof_indices (local_dof_indices); - hanging_node_constraints.distribute_local_to_global (cell_matrix, - local_dof_indices, - system_matrix); - } - - timer.leave_subsection(); -} - - - -template -void HelmholtzProblem::solve () -{ - for (unsigned int i=0; i -void HelmholtzProblem::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 fe(element_degree); - HelmholtzProblem problem(fe); - problem.run(); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} - - - diff --git a/tests/benchmarks/testlist.sh b/tests/benchmarks/testlist.sh deleted file mode 100755 index 96858aeac9..0000000000 --- a/tests/benchmarks/testlist.sh +++ /dev/null @@ -1,2 +0,0 @@ -#!/bin/bash -export TESTS="step-22 tablehandler test_assembly test_poisson test_hp" -- 2.39.5