From: oliver Date: Tue, 20 Dec 2005 20:12:07 +0000 (+0000) Subject: Added a first version of an example, which shows the use of the new X-Git-Tag: v8.0.0~12757 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c62189d0e689466a45137b9ab419039787b9ac2a;p=dealii.git Added a first version of an example, which shows the use of the new hpDoFHandler functionality. git-svn-id: https://svn.dealii.org/trunk@11898 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-21/Makefile b/deal.II/examples/step-21/Makefile new file mode 100644 index 0000000000..5cecab44d4 --- /dev/null +++ b/deal.II/examples/step-21/Makefile @@ -0,0 +1,154 @@ +# $Id$ + + +# For the small projects Makefile, you basically need to fill in only +# four fields. +# +# The first is the name of the application. It is assumed that the +# application name is the same as the base file name of the single C++ +# file from which the application is generated. +target = $(basename $(shell echo step-*.cc)) + +# The second field determines whether you want to run your program in +# debug or optimized mode. The latter is significantly faster, but no +# run-time checking of parameters and internal states is performed, so +# you should set this value to `on' while you develop your program, +# and to `off' when running production computations. +debug-mode = on + + +# As third field, we need to give the path to the top-level deal.II +# directory. You need to adjust this to your needs. Since this path is +# probably the most often needed one in the Makefile internals, it is +# designated by a single-character variable, since that can be +# reference using $D only, i.e. without the parentheses that are +# required for most other parameters, as e.g. in $(target). +D = ../../ + + +# The last field specifies the names of data and other files that +# shall be deleted when calling `make clean'. Object and backup files, +# executables and the like are removed anyway. Here, we give a list of +# files in the various output formats that deal.II supports. +clean-up-files = *gmv *gnuplot *gpl *eps *pov + + + + +# +# +# Usually, you will not need to change something beyond this point. +# +# +# The next statement tell the `make' program where to find the +# deal.II top level directory and to include the file with the global +# settings +include $D/common/Make.global_options + + +# Since the whole project consists of only one file, we need not +# consider difficult dependencies. We only have to declare the +# libraries which we want to link to the object file, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-2d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# We now use the variable defined above which switch between debug and +# optimized mode to select the set of libraries to link with. Included +# in the list of libraries is the name of the object file which we +# will produce from the single C++ file. Note that by default we use +# the extension .g.o for object files compiled in debug mode and .o for +# object files in optimized mode (or whatever the local default on your +# system is instead of .o). +ifeq ($(debug-mode),on) + libraries = $(target).g.$(OBJEXT) $(libs.g) +else + libraries = $(target).$(OBJEXT) $(libs.o) +endif + + +# Now comes the first production rule: how to link the single object +# file produced from the single C++ file into the executable. Since +# this is the first rule in the Makefile, it is the one `make' selects +# if you call it without arguments. +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) + + +# To make running the application somewhat independent of the actual +# program name, we usually declare a rule `run' which simply runs the +# program. You can then run it by typing `make run'. This is also +# useful if you want to call the executable with arguments which do +# not change frequently. You may then want to add them to the +# following rule: +run: $(target) + @echo ============================ Running $< + @./$(target)$(EXEEXT) + + +# As a last rule to the `make' program, we define what to do when +# cleaning up a directory. This usually involves deleting object files +# and other automatically created files such as the executable itself, +# backup files, and data files. Since the latter are not usually quite +# diverse, you needed to declare them at the top of this file. +clean: + -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) + + +# Since we have not yet stated how to make an object file from a C++ +# file, we should do so now. Since the many flags passed to the +# compiler are usually not of much interest, we suppress the actual +# command line using the `at' sign in the first column of the rules +# and write the string indicating what we do instead. +./%.g.$(OBJEXT) : + @echo ==============debug========= $( $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then : else rm $@ ; fi + + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc new file mode 100644 index 0000000000..3223c43b58 --- /dev/null +++ b/deal.II/examples/step-21/step-21.cc @@ -0,0 +1,1825 @@ +/* $Id$ */ +/* Authors: Ralf Hartmann, University of Heidelberg, 2001 */ +/* hp-Version, Oliver Kayser-Herold, TU-Braunschweig 2005 */ + +/* $Id$ */ +/* Version: $Name$ */ +/* */ +/* Copyright (C) 2001-2005 by the deal.II authors */ +/* */ +/* This file is subject to QPL and may not be distributed */ +/* without copyright and license information. Please refer */ +/* to the file deal.II/doc/license.html for the text and */ +/* further information on this license. */ + + // 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 + // Instead of the usual DoFHandler + // the hpDoFHandler class has to be + // used to gain access to the + // hp-Functionality. The hpDoFHandler + // provides essentially the same + // interface as the standard DoFHandler. +#include + // Due to the implementation of the + // hp-Method, the DoFAccessor classes + // stay the same as for the DoFHandler, + // and hence also require the same + // include files. +#include +#include +#include + +#include +#include + + // The PETSc vector and matrix + // classes are used to have access to + // the PETSc preconditioners. +#include +#include + // Then we also need interfaces for solvers + // and preconditioners that PETSc provides: +#include +#include + + // We are going to use gradients as + // refinement indicator. +#include + // Finally we do some time comparison + // using the ``Timer'' class. +#include + + // And this again is C++: +#include +#include + + // Now some additional utility classes will + // be necessary. As the name "collection" + // suggests, these are essentially + // container classes which store the + // quadrature and finite element objects, + // for the different polynomial degrees. +#include +#include + + // The first of the following + // two files provides the hpFEValues + // class, which implements the same + // functionality as the FEValues class + // with the difference that it takes + // "collection" objects instead of + // single finite element or quadrature + // objects. I.e. instead of the + // usual Quadrature object a + // QCollection object is needed. The + // select_fe_values file finally + // provides some wrappers, which should + // enable the user to write DoFHandler + // independent code. +#include +#include + + // A compressed sparsity pattern is + // not an explicit prerequisite for the + // use of the hp-functionality. But as + // a standard sparsity pattern has to + // based on a bandwidth estimate for the + // highest polynomial degree, it would + // be prohibitively bad. Therefore, + // the recommended way is to explicitly + // build a compressed sparsity pattern before + // creating the matrices. +#include + + + // @sect3{Equation data} + // + // First we define the classes + // representing the equation-specific + // functions. Both classes, ``RHS'' + // and ``BoundaryValues'', are + // derived from the ``Function'' + // class. Only the ``value_list'' + // function are implemented because + // only lists of function values are + // computed rather than single + // values. +template +class RHS: public Function +{ + public: + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component=0) const; +}; + + +template +class BoundaryValues: public Function +{ + public: + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component=0) const; +}; + + + // The class ``Beta'' represents the + // vector valued flow field of the + // linear transport equation and is + // not derived from the ``Function'' + // class as we prefer to get function + // values of type ``Point'' rather + // than of type + // ``Vector''. This, because + // there exist scalar products + // between ``Point'' and ``Point'' as + // well as between ``Point'' and + // ``Tensor'', simplifying terms like + // $\beta\cdot n$ and + // $\beta\cdot\nabla v$. + // + // An unnecessary empty constructor + // is added to the class to work + // around a bug in Compaq's cxx + // compiler which otherwise reports + // an error about an omitted + // initializer for an object of + // this class further down. +template +class Beta +{ + public: + Beta () {}; + void value_list (const std::vector > &points, + std::vector > &values) const; +}; + + + // The implementation of the + // ``value_list'' functions of these + // classes are rather simple. For + // simplicity the right hand side is + // set to be zero but will be + // assembled anyway. +template +void RHS::value_list(const std::vector > &points, + std::vector &values, + const unsigned int) const +{ + // Usually we check whether input + // parameter have the right sizes. + Assert(values.size()==points.size(), + ExcDimensionMismatch(values.size(),points.size())); + + for (unsigned int i=0; i +void Beta::value_list(const std::vector > &points, + std::vector > &values) const +{ + Assert(values.size()==points.size(), + ExcDimensionMismatch(values.size(),points.size())); + + for (unsigned int i=0; i &p=points[i]; + Point &beta=values[i]; + + beta(0) = -p(1); + beta(1) = p(0); + beta /= std::sqrt(beta.square()); + } +} + + + // Hence the inflow boundary of the + // unit square [0,1]^2 are the right + // and the lower boundaries. We + // prescribe discontinuous boundary + // values 1 and 0 on the x-axis and + // value 0 on the right boundary. The + // values of this function on the + // outflow boundaries will not be + // used within the DG scheme. +template +void BoundaryValues::value_list(const std::vector > &points, + std::vector &values, + const unsigned int) const +{ + Assert(values.size()==points.size(), + ExcDimensionMismatch(values.size(),points.size())); + + for (unsigned int i=0; i +class DGTransportEquation +{ + public: + DGTransportEquation(); + + void assemble_cell_term(const FEValues& fe_v, + FullMatrix &u_v_matrix, + Vector &cell_vector) const; + + void assemble_boundary_term(const FEFaceValues& fe_v, + FullMatrix &u_v_matrix, + Vector &cell_vector) const; + + void assemble_face_term1(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix) const; + + void assemble_face_term2(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix, + FullMatrix &u_vn_matrix, + FullMatrix &un_vn_matrix) const; + private: + const Beta beta_function; + const RHS rhs_function; + const BoundaryValues boundary_function; +}; + + +template +DGTransportEquation::DGTransportEquation () + : + beta_function (), + rhs_function (), + boundary_function () +{} + + + // @sect4{Function: assemble_cell_term} + // + // The ``assemble_cell_term'' + // function assembles the cell terms + // of the discretization. + // ``u_v_matrix'' is a cell matrix, + // i.e. for a DG method of degree 1, + // it is of size 4 times 4, and + // ``cell_vector'' is of size 4. + // When this function is invoked, + // ``fe_v'' is already reinit'ed with the + // current cell before and includes + // all shape values needed. +template +void DGTransportEquation::assemble_cell_term( + const FEValues &fe_v, + FullMatrix &u_v_matrix, + Vector &cell_vector) const +{ + // First we ask ``fe_v'' for the + // shape gradients, shape values and + // quadrature weights, + const std::vector &JxW = fe_v.get_JxW_values (); + + // Then the flow field beta and the + // ``rhs_function'' are evaluated at + // the quadrature points, + std::vector > beta (fe_v.n_quadrature_points); + std::vector rhs (fe_v.n_quadrature_points); + + beta_function.value_list (fe_v.get_quadrature_points(), beta); + rhs_function.value_list (fe_v.get_quadrature_points(), rhs); + + // and the cell matrix and cell + // vector are assembled due to the + // terms $-(u,\beta\cdot\nabla + // v)_K$ and $(f,v)_K$. + for (unsigned int point=0; point +void DGTransportEquation::assemble_boundary_term( + const FEFaceValues& fe_v, + FullMatrix &u_v_matrix, + Vector &cell_vector) const +{ + // Again, as in the previous + // function, we ask the ``FEValues'' + // object for the shape values and + // the quadrature weights + const std::vector &JxW = fe_v.get_JxW_values (); + // but here also for the normals. + const std::vector > &normals = fe_v.get_normal_vectors (); + + // We evaluate the flow field + // and the boundary values at the + // quadrature points. + std::vector > beta (fe_v.n_quadrature_points); + std::vector g(fe_v.n_quadrature_points); + + beta_function.value_list (fe_v.get_quadrature_points(), beta); + boundary_function.value_list (fe_v.get_quadrature_points(), g); + + // Then we assemble cell vector and + // cell matrix according to the DG + // method given in the + // introduction. + for (unsigned int point=0; point0) + for (unsigned int i=0; i +void DGTransportEquation::assemble_face_term1( + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix) const +{ + // Again, as in the previous + // function, we ask the FEValues + // objects for the shape values, + // the quadrature weights and the + // normals + const std::vector &JxW = fe_v.get_JxW_values (); + const std::vector > &normals = fe_v.get_normal_vectors (); + + // and we evaluate the flow field + // at the quadrature points. + std::vector > beta (fe_v.n_quadrature_points); + beta_function.value_list (fe_v.get_quadrature_points(), beta); + + // Then we assemble the cell + // matrices according to the DG + // method given in the + // introduction. + for (unsigned int point=0; point0) + for (unsigned int i=0; i +void DGTransportEquation::assemble_face_term2( + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix, + FullMatrix &u_vn_matrix, + FullMatrix &un_vn_matrix) const +{ + // the first few lines are the same + const std::vector &JxW = fe_v.get_JxW_values (); + const std::vector > &normals = fe_v.get_normal_vectors (); + + std::vector > beta (fe_v.n_quadrature_points); + + beta_function.value_list (fe_v.get_quadrature_points(), beta); + + for (unsigned int point=0; point0) + { + // This terms we've already seen. + for (unsigned int i=0; i +class DGMethod +{ + public: + DGMethod (); + ~DGMethod (); + + void run (); + + private: + void setup_system (); + void assemble_system1 (); + void assemble_system2 (); + void solve (PETScWrappers::Vector &solution); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + + // Currently this example code + // works on a rectangular domain. + // Hence a linear mapping of degree + // 1 is sufficient. If irregular + // boundaries come into play, this + // object should be replaced by + // a MappingCollection. + const MappingQ1 mapping; + + // In contrast to the example code + // of step-12, this time DG elements + // of different degree will be used. + // The different FiniteElement + // objects for the different polynomial + // degrees will be stored in the + // fe_collection object. + FECollection fe_collection; + + // As already mentioned, the + // standard DoFHandler has to be + // replaced by a hpDoFHandler. + hpDoFHandler dof_handler; + + PETScWrappers::SparseMatrix system_matrix; + + // We define the quadrature + // formulae for the cell and the + // face terms of the + // discretization. + // Clearly the hp-Method requires + // a complete set of quadrature + // rules for each polynomial + // degree which will be used in the + // computations. + std::vector *> quad; + QCollection *quadrature; + std::vector *> face_quad; + QCollection *face_quadrature; + + // And there are two solution + // vectors, that store the + // solutions to the problems + // corresponding to the two + // different assembling routines + // ``assemble_system1'' and + // ``assemble_system2''; + PETScWrappers::Vector solution1; + PETScWrappers::Vector solution2; + PETScWrappers::Vector right_hand_side; + + // Finally this class includes an + // object of the + // DGTransportEquations class + // described above. + const DGTransportEquation dg; +}; + + +template +DGMethod::DGMethod () + : + mapping (), + dof_handler (triangulation), + dg () +{ + quadrature = new QCollection (); + face_quadrature = new QCollection (); + + Quadrature *quad_temp; + Quadrature *face_quad_temp; + // Change here for hp + // methods of + // different maximum degrees. + const unsigned int hp_degree = 5; + unsigned int elm_no; + for (unsigned int i = 0; i < hp_degree; ++i) + { + elm_no = fe_collection.add_fe (FE_DGQ (i+1)); + + quad_temp = new QGauss (i+2); + quad.push_back (quad_temp); + quadrature->add_quadrature (*quad_temp); + + //TODO: This is currently a design error in the hp Framework. + // While the idea of using a QuadratureCollection for the + // element stiffness matrices is straightforward, it is + // not complete for the face terms. At these places, we can + // have different polynomial degrees on the two faces. Therefore + // a quadrature rule, which is sufficient for the highest polynomial + // degree on both faces, would be required. This will certainly need + // some changes in the hpFEValues classes. But for a first example, + // simply using a uniformly high quadrature rule will do the job + // although is obviously not very efficient, as even linear polynomials + // will be integrated with 7th order Gauss formulas. + face_quad_temp = new QGauss (7); // (i+2); + face_quad.push_back (face_quad_temp); + face_quadrature->add_quadrature (*face_quad_temp); + } +} + + +template +DGMethod::~DGMethod () +{ + dof_handler.clear (); + + // Free the dynamically created quadrature rules + delete quadrature; + delete face_quadrature; + + typename std::vector *>::iterator quad_iter; + for (quad_iter = quad.begin (); quad_iter != quad.end (); ++quad_iter) + delete (*quad_iter); + + typename std::vector *>::iterator face_quad_iter; + for (face_quad_iter = face_quad.begin (); + face_quad_iter != face_quad.end (); ++face_quad_iter) + delete (*face_quad_iter); +} + + +template +void DGMethod::setup_system () +{ + // First we need to distribute the + // DoFs. + dof_handler.distribute_dofs (fe_collection); + + // The DoFs of a cell are coupled + // with all DoFs of all neighboring + // cells. Therefore the maximum + // number of matrix entries per row + // is needed when all neighbors of + // a cell are once more refined + // than the cell under + // consideration. + CompressedSparsityPattern compressed_pattern (dof_handler.n_dofs ()); + DoFTools::make_sparsity_pattern (dof_handler, compressed_pattern); + DoFTools::make_flux_sparsity_pattern (dof_handler, compressed_pattern); + + system_matrix.reinit (compressed_pattern); + + solution1.reinit (dof_handler.n_dofs()); + solution2.reinit (dof_handler.n_dofs()); + right_hand_side.reinit (dof_handler.n_dofs()); +} + + + // @sect4{Function: assemble_system1} + // + // We proceed with the + // ``assemble_system1'' function that + // implements the DG discretization + // in its first version. This + // function repeatedly calls the + // ``assemble_cell_term'', + // ``assemble_boundary_term'' and + // ``assemble_face_term1'' functions + // of the ``DGTransportEquation'' + // object. The + // ``assemble_boundary_term'' covers + // the first case mentioned in the + // introduction. + // + // 1. face is at boundary + // + // This function takes a + // ``FEFaceValues'' object as + // argument. In contrast to that + // ``assemble_face_term1'' + // takes two ``FEFaceValuesBase'' + // objects; one for the shape + // functions on the current cell and + // the other for shape functions on + // the neighboring cell under + // consideration. Both objects are + // either of class ``FEFaceValues'' + // or of class ``FESubfaceValues'' + // (both derived from + // ``FEFaceValuesBase'') according to + // the remaining cases mentioned + // in the introduction: + // + // 2. neighboring cell is finer + // (current cell: ``FESubfaceValues'', + // neighboring cell: ``FEFaceValues''); + // + // 3. neighboring cell is of the same + // refinement level (both, current + // and neighboring cell: + // ``FEFaceValues''); + // + // 4. neighboring cell is coarser + // (current cell: ``FEFaceValues'', + // neighboring cell: + // ``FESubfaceValues''). + // + // If we considered globally refined + // meshes then only case 3 would + // occur. But as we consider also + // locally refined meshes we need to + // distinguish all four cases making + // the following assembling function + // a bit longish. +template +void DGMethod::assemble_system1 () +{ + // First we create the + // ``UpdateFlags'' for the + // ``FEValues'' and the + // ``FEFaceValues'' objects. + UpdateFlags update_flags = update_values + | update_gradients + | update_q_points + | update_JxW_values; + + // Note, that on faces we do not + // need gradients but we need + // normal vectors. + UpdateFlags face_update_flags = update_values + | update_q_points + | update_JxW_values + | update_normal_vectors; + + // On the neighboring cell we only + // need the shape values. Given a + // specific face, the quadrature + // points and `JxW values' are the + // same as for the current cells, + // the normal vectors are known to + // be the negative of the normal + // vectors of the current cell. + UpdateFlags neighbor_face_update_flags = update_values; + + // Then we create the ``FEValues'' + // object. Note, that since version + // 3.2.0 of deal.II the constructor + // of this class takes a + // ``Mapping'' object as first + // argument. Although the + // constructor without ``Mapping'' + // argument is still supported it + // is recommended to use the new + // constructor. This reduces the + // effect of `hidden magic' (the + // old constructor implicitely + // assumes a ``MappingQ1'' mapping) + // and makes it easier to change + // the mapping object later. + typename SelectFEValues >::FEValues fe_v_x ( + mapping, fe_collection, *quadrature, update_flags); + + // Similarly we create the + // ``FEFaceValues'' and + // ``FESubfaceValues'' objects for + // both, the current and the + // neighboring cell. Within the + // following nested loop over all + // cells and all faces of the cell + // they will be reinited to the + // current cell and the face (and + // subface) number. + typename SelectFEValues >::FEFaceValues fe_v_face_x ( + mapping, fe_collection, *face_quadrature, face_update_flags); + typename SelectFEValues >::FESubfaceValues fe_v_subface_x ( + mapping, fe_collection, *face_quadrature, face_update_flags); + typename SelectFEValues >::FEFaceValues fe_v_face_neighbor_x ( + mapping, fe_collection, *face_quadrature, neighbor_face_update_flags); + typename SelectFEValues >::FESubfaceValues fe_v_subface_neighbor_x ( + mapping, fe_collection, *face_quadrature, neighbor_face_update_flags); + + // Now we create the cell matrices + // and vectors. Here we need two + // cell matrices, both for face + // terms that include test + // functions ``v'' (shape functions + // of the current cell). To be more + // precise, the first matrix will + // include the `u and v terms' and + // the second that will include the + // `un and v terms'. Here we recall + // the convention that `un' is + // the shortcut for `u_neighbor' + // and represents the $u_hat$, see + // introduction. + const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell (); + + FullMatrix u_v_matrix (max_dofs_per_cell, max_dofs_per_cell); + FullMatrix un_v_matrix (max_dofs_per_cell, max_dofs_per_cell); + Vector cell_vector (max_dofs_per_cell); + + // Furthermore we need some cell + // iterators. + typename hpDoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + // Now we start the loop over all + // active cells. + for (;cell!=endc; ++cell) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + std::vector dofs_neighbor; + + // In the + // ``assemble_face_term1'' + // function contributions to + // the cell matrices and the + // cell vector are only + // ADDED. Therefore on each + // cell we need to reset the + // ``u_v_matrix'' and + // ``cell_vector'' to zero, + // before assembling the cell terms. + u_v_matrix = 0; + cell_vector = 0; + + // Now we reinit the ``FEValues'' + // object for the current cell + fe_v_x.reinit (cell); + + // and call the function + // that assembles the cell + // terms. The first argument is + // the ``FEValues'' that was + // previously reinit'ed on the + // current cell. + dg.assemble_cell_term(fe_v_x.get_present_fe_values (), + u_v_matrix, + cell_vector); + + // As in previous examples the + // vector `dofs' includes the + // dof_indices. + dofs.resize (dofs_per_cell); + cell->get_dof_indices (dofs); + + // This is the start of the + // nested loop over all faces. + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + // First we set the face + // iterator + typename hpDoFHandler::face_iterator face=cell->face(face_no); + + // and clear the + // ``un_v_matrix'' on each + // face. + un_v_matrix = 0; + + // Now we distinguish the + // four different cases in + // the ordering mentioned + // above. We start with + // faces belonging to the + // boundary of the domain. + if (face->at_boundary()) + { + // We reinit the + // ``FEFaceValues'' + // object to the + // current face + fe_v_face_x.reinit (cell, face_no); + + // and assemble the + // corresponding face + // terms. + dg.assemble_boundary_term(fe_v_face_x.get_present_fe_values (), + u_v_matrix, + cell_vector); + } + else + { + // Now we are not on + // the boundary of the + // domain, therefore + // there must exist a + // neighboring cell. + typename hpDoFHandler::cell_iterator neighbor= + cell->neighbor(face_no); + + // We proceed with the + // second and most + // complicated case: + // the neighboring cell + // is more refined than + // the current cell. As + // in deal.II + // neighboring cells + // are restricted to + // have a level + // difference of not + // more than one, the + // neighboring cell is + // known to be at most + // ONCE more refined + // than the current + // cell. Furthermore + // also the face is + // more refined, + // i.e. it has + // children. Here we + // note that the + // following part of + // code will not work + // for ``dim==1''. + if (face->has_children()) + { + // First we store + // which number the + // current cell has + // in the list of + // neighbors of the + // neighboring + // cell. Hence, + // neighbor->neighbor(neighbor2) + // equals the + // current cell + // ``cell''. + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face_no); + + + // We loop over + // subfaces + for (unsigned int subface_no=0; + subface_no::subfaces_per_face; + ++subface_no) + { + // and set the + // cell + // iterator + // ``neighbor_child'' + // to the cell + // placed + // `behind' the + // current + // subface. + typename hpDoFHandler::active_cell_iterator neighbor_child= + neighbor->child(GeometryInfo:: + child_cell_on_face(neighbor2,subface_no)); + + // As these are + // quite + // complicated + // indirections + // which one + // does not + // usually get + // right at + // first + // attempt we + // check for + // the internal + // consistency. + Assert (neighbor_child->face(neighbor2) == face->child(subface_no), + ExcInternalError()); + Assert (!neighbor_child->has_children(), ExcInternalError()); + + // We need to + // reset the + // ``un_v_matrix'' + // on each + // subface + // because on + // each subface + // the ``un'' + // belong to + // different + // neighboring + // cells. + un_v_matrix = 0; + + // As already + // mentioned + // above for + // the current + // case (case + // 2) we employ + // the + // ``FESubfaceValues'' + // of the + // current + // cell (here + // reinited for + // the current + // cell, face + // and subface) + // and we + // employ the + // FEFaceValues + // of the + // neighboring + // child cell. + fe_v_subface_x.reinit (cell, face_no, subface_no); + fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2); + + dg.assemble_face_term1(fe_v_subface_x.get_present_fe_values (), + fe_v_face_neighbor_x.get_present_fe_values (), + u_v_matrix, + un_v_matrix); + + // Then we get + // the dof + // indices of + // the + // neighbor_child + // cell + dofs_neighbor.resize (neighbor_child->get_fe().dofs_per_cell); + neighbor_child->get_dof_indices (dofs_neighbor); + + // and + // distribute + // ``un_v_matrix'' + // to the + // system_matrix + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + for (unsigned int k=0; kget_fe().dofs_per_cell; ++k) + system_matrix.add(dofs[i], dofs_neighbor[k], + un_v_matrix(i,k)); + } + // End of ``if + // (face->has_children())'' + } + else + { + // We proceed with + // case 3, + // i.e. neighboring + // cell is of the + // same refinement + // level as the + // current cell. + if (neighbor->level() == cell->level()) + { + // Like before + // we store + // which number + // the current + // cell has in + // the list of + // neighbors of + // the + // neighboring + // cell. + const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); + + // We reinit + // the + // ``FEFaceValues'' + // of the + // current and + // neighboring + // cell to the + // current face + // and assemble + // the + // corresponding + // face terms. + fe_v_face_x.reinit (cell, face_no); + fe_v_face_neighbor_x.reinit (neighbor, neighbor2); + + dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (), + fe_v_face_neighbor_x.get_present_fe_values (), + u_v_matrix, + un_v_matrix); + // End of ``if + // (neighbor->level() + // == + // cell->level())'' + } + else + { + // Finally we + // consider + // case 4. When + // the + // neighboring + // cell is not + // finer and + // not of the + // same + // refinement + // level as the + // current cell + // it must be + // coarser. + Assert(neighbor->level() < cell->level(), ExcInternalError()); + + // Find out the + // how many'th + // face_no and + // subface_no + // the current + // face is + // w.r.t. the + // neighboring + // cell. + const std::pair faceno_subfaceno= + cell->neighbor_of_coarser_neighbor(face_no); + const unsigned int neighbor_face_no=faceno_subfaceno.first, + neighbor_subface_no=faceno_subfaceno.second; + + Assert (neighbor->neighbor(neighbor_face_no) + ->child(GeometryInfo::child_cell_on_face( + face_no,neighbor_subface_no)) == cell, ExcInternalError()); + + // Reinit the + // appropriate + // ``FEFaceValues'' + // and assemble + // the face + // terms. + fe_v_face_x.reinit (cell, face_no); + fe_v_subface_neighbor_x.reinit (neighbor, neighbor_face_no, + neighbor_subface_no); + + dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (), + fe_v_subface_neighbor_x.get_present_fe_values (), + u_v_matrix, + un_v_matrix); + } + + // Now we get the + // dof indices of + // the + // ``neighbor_child'' + // cell, + dofs_neighbor.resize (neighbor->get_fe().dofs_per_cell); + neighbor->get_dof_indices (dofs_neighbor); + + // and distribute the + // ``un_v_matrix''. + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + for (unsigned int k=0; kget_fe().dofs_per_cell; ++k) + system_matrix.add(dofs[i], dofs_neighbor[k], + un_v_matrix(i,k)); + } + // End of ``face not at boundary'': + } + // End of loop over all faces: + } + + // Finally we distribute the + // ``u_v_matrix'' + for (unsigned int i=0; i +void DGMethod::assemble_system2 () +{ + UpdateFlags update_flags = update_values + | update_gradients + | update_q_points + | update_JxW_values; + + UpdateFlags face_update_flags = update_values + | update_q_points + | update_JxW_values + | update_normal_vectors; + + UpdateFlags neighbor_face_update_flags = update_values; + + // Here we do not need + // ``fe_v_face_neighbor'' as case 4 + // does not occur. + typename SelectFEValues >::FEValues fe_v_x ( + mapping, fe_collection, *quadrature, update_flags); + typename SelectFEValues >::FEFaceValues fe_v_face_x ( + mapping, fe_collection, *face_quadrature, face_update_flags); + typename SelectFEValues >::FESubfaceValues fe_v_subface_x ( + mapping, fe_collection, *face_quadrature, face_update_flags); + typename SelectFEValues >::FEFaceValues fe_v_face_neighbor_x ( + mapping, fe_collection, *face_quadrature, neighbor_face_update_flags); + + const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell (); + + FullMatrix u_v_matrix (max_dofs_per_cell, max_dofs_per_cell); + FullMatrix un_v_matrix (max_dofs_per_cell, max_dofs_per_cell); + + // Additionally we need the + // following two cell matrices, + // both for face term that include + // test function ``vn'' (shape + // functions of the neighboring + // cell). To be more precise, the + // first matrix will include the `u + // and vn terms' and the second + // that will include the `un and vn + // terms'. + FullMatrix u_vn_matrix (max_dofs_per_cell, max_dofs_per_cell); + FullMatrix un_vn_matrix (max_dofs_per_cell, max_dofs_per_cell); + + Vector cell_vector (max_dofs_per_cell); + + // The following lines are roughly + // the same as in the previous + // function. + typename hpDoFHandler::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; + std::vector dofs (dofs_per_cell); + std::vector dofs_neighbor (dofs_per_cell); + + u_v_matrix = 0; + cell_vector = 0; + + fe_v_x.reinit (cell); + + dg.assemble_cell_term(fe_v_x.get_present_fe_values (), + u_v_matrix, + cell_vector); + + cell->get_dof_indices (dofs); + + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + typename hpDoFHandler::face_iterator face= + cell->face(face_no); + + // Case 1: + if (face->at_boundary()) + { + fe_v_face_x.reinit (cell, face_no); + + dg.assemble_boundary_term(fe_v_face_x.get_present_fe_values (), + u_v_matrix, + cell_vector); + } + else + { + Assert (cell->neighbor(face_no).state() == IteratorState::valid, + ExcInternalError()); + typename hpDoFHandler::cell_iterator neighbor= + cell->neighbor(face_no); + + const unsigned int dofs_on_neighbor = neighbor->get_fe().dofs_per_cell; + + // Case 2: + if (face->has_children()) + { + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face_no); + + for (unsigned int subface_no=0; + subface_no::subfaces_per_face; + ++subface_no) + { + typename hpDoFHandler::cell_iterator neighbor_child= + neighbor->child(GeometryInfo::child_cell_on_face( + neighbor2,subface_no)); + const unsigned int dofs_on_neighbor_child = neighbor_child->get_fe().dofs_per_cell; + + Assert (neighbor_child->face(neighbor2) == face->child(subface_no), + ExcInternalError()); + Assert (!neighbor_child->has_children(), ExcInternalError()); + + un_v_matrix = 0; + u_vn_matrix = 0; + un_vn_matrix = 0; + + fe_v_subface_x.reinit (cell, face_no, subface_no); + fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2); + + dg.assemble_face_term2(fe_v_subface_x.get_present_fe_values (), + fe_v_face_neighbor_x.get_present_fe_values (), + u_v_matrix, + un_v_matrix, + u_vn_matrix, + un_vn_matrix); + + dofs_neighbor.resize (dofs_on_neighbor_child); + neighbor_child->get_dof_indices (dofs_neighbor); + + for (unsigned int i=0; ilevel() == cell->level() && + neighbor->index() > cell->index()) + { + const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); + + un_v_matrix = 0; + u_vn_matrix = 0; + un_vn_matrix = 0; + + fe_v_face_x.reinit (cell, face_no); + fe_v_face_neighbor_x.reinit (neighbor, neighbor2); + + dg.assemble_face_term2(fe_v_face_x.get_present_fe_values (), + fe_v_face_neighbor_x.get_present_fe_values (), + u_v_matrix, + un_v_matrix, + u_vn_matrix, + un_vn_matrix); + + dofs_neighbor.resize (dofs_on_neighbor); + neighbor->get_dof_indices (dofs_neighbor); + + for (unsigned int i=0; i +void DGMethod::solve (PETScWrappers::Vector &solution) +{ + SolverControl solver_control (solution.size(), + 1e-13, false, false); + PETScWrappers::SolverBiCG bicg (solver_control); + PETScWrappers::PreconditionILU preconditioner(system_matrix); + + // Then solve the system: + bicg.solve (system_matrix, solution, right_hand_side, + preconditioner); + + solution.compress (); + + std::cout << "Iterations : " << solver_control.last_step() << std::endl; +} + + + // We refine the grid according to a + // very simple refinement criterion, + // namely an approximation to the + // gradient of the solution. As here + // we consider the DG(1) method + // (i.e. we use piecewise bilinear + // shape functions) we could simply + // compute the gradients on each + // cell. But we do not want to base + // our refinement indicator on the + // gradients on each cell only, but + // want to base them also on jumps of + // the discontinuous solution + // function over faces between + // neighboring cells. The simpliest + // way of doing that is to compute + // approximative gradients by + // difference quotients including the + // cell under consideration and its + // neighbors. This is done by the + // ``DerivativeApproximation'' class + // that computes the approximate + // gradients in a way similar to the + // ``GradientEstimation'' described + // in Step 9 of this tutorial. In + // fact, the + // ``DerivativeApproximation'' class + // was developed following the + // ``GradientEstimation'' class of + // Step 9. Relating to the + // discussion in Step 9, here we + // consider $h^{1+d/2}|\nabla_h + // u_h|$. Futhermore we note that we + // do not consider approximate second + // derivatives because solutions to + // the linear advection equation are + // in general not in H^2 but in H^1 + // (to be more precise, in H^1_\beta) + // only. +template +void DGMethod::refine_grid () +{ + // The ``DerivativeApproximation'' + // class computes the gradients to + // float precision. This is + // sufficient as they are + // approximate and serve as + // refinement indicators only. + Vector gradient_indicator (triangulation.n_active_cells()); + + // Now the approximate gradients + // are computed + DerivativeApproximation::approximate_gradient (mapping, + dof_handler, + solution2, + gradient_indicator); + + // and they are cell-wise scaled by + // the factor $h^{1+d/2}$ + typename hpDoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + gradient_indicator(cell_no)*=std::pow(cell->diameter(), 1+1.0*dim/2); + + // Finally they serve as refinement + // indicator. + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + gradient_indicator, + 0.3, 0.1); + + // Simple heuristic hp-Refinement. + // As the indicator marks nonsmooth + // regions, p-refine all non marked + // regions, while the marked + // regions clearly deserve an + // h-refinement. + cell = dof_handler.begin_active (); + for (; cell!=endc; ++cell) + if (!cell->refine_flag_set ()) + cell->set_active_fe_index (cell->active_fe_index () + 1); + + triangulation.execute_coarsening_and_refinement (); +} + + + // The output of this program + // consists of eps-files of the + // adaptively refined grids and the + // numerical solutions given in + // gnuplot format. This was covered + // in previous examples and will not + // be further commented on. +template +void DGMethod::output_results (const unsigned int cycle) const +{ + // Write the grid in eps format. + std::string filename = "grid-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".eps"; + std::cout << "Writing grid to <" << filename << ">..." << std::endl; + std::ofstream eps_output (filename.c_str()); + + GridOut grid_out; + grid_out.write_eps (triangulation, eps_output); + + // Output of the solution in + // gnuplot format. + filename = "sol-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + // filename += ".gnuplot"; + filename += ".vtk"; + std::cout << "Writing solution to <" << filename << ">..." + << std::endl << std::endl; + std::ofstream gnuplot_output (filename.c_str()); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution2, "u"); + + data_out.build_patches (5); + +// data_out.write_gnuplot(gnuplot_output); + data_out.write_vtk(gnuplot_output); +} + + + // The following ``run'' function is + // similar to previous examples. The + // only difference is that the + // problem is assembled and solved + // twice on each refinement step; + // first by ``assemble_system1'' that + // implements the first version and + // then by ``assemble_system2'' that + // implements the second version of + // writing the DG + // discretization. Furthermore the + // time needed by each of the two + // assembling routines is measured. +template +void DGMethod::run () +{ + for (unsigned int cycle=0; cycle<5; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + + triangulation.refine_global (3); + } + else + refine_grid (); + + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + // The constructor of the Timer + // class automatically starts + // the time measurement. + Timer assemble_timer; + // First assembling routine. + assemble_system1 (); + // The operator () accesses the + // current time without + // disturbing the time + // measurement. + std::cout << "Time of assemble_system1: " + << assemble_timer() + << std::endl; + solve (solution1); + + + // As preparation for the + // second assembling routine we + // reinit the system matrix, the + // right hand side vector and + // the Timer object. + system_matrix = 0; + right_hand_side = 0; + assemble_timer.reset(); + + // We start the Timer, + assemble_timer.start(); + // call the second assembling routine + assemble_system2 (); + // and access the current time. + std::cout << "Time of assemble_system2: " + << assemble_timer() + << std::endl; + solve (solution2); + + // To make sure that both + // versions of the DG method + // yield the same + // discretization and hence the + // same solution we check the + // two solutions for equality. + solution1-=solution2; + + const double difference=solution1.linfty_norm(); + if (difference>1e-13) + std::cout << "solution1 and solution2 differ!!" << std::endl; + else + std::cout << "solution1 and solution2 coincide." << std::endl; + + // Finally we perform the + // output. + output_results (cycle); + } +} + + // The following ``main'' function is + // similar to previous examples and + // need not to be commented on. +int main (int argc, char **argv) +{ + try + { + PetscInitialize(&argc,&argv,0,0); + + { + DGMethod<2> dgmethod; + dgmethod.run (); + } + + PetscFinalize (); + } + 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; +} + +