From: Guido Kanschat Date: Fri, 18 Sep 2009 18:38:23 +0000 (+0000) Subject: add example 38, based on 12 X-Git-Tag: v8.0.0~7100 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=41d926e8797e963573af7281502453accd23fd38;p=dealii.git add example 38, based on 12 git-svn-id: https://svn.dealii.org/trunk@19475 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-38/Makefile b/deal.II/examples/step-38/Makefile new file mode 100644 index 0000000000..db1de21066 --- /dev/null +++ b/deal.II/examples/step-38/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 anything 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-38/doc/intro.dox b/deal.II/examples/step-38/doc/intro.dox new file mode 100644 index 0000000000..2f9248cc02 --- /dev/null +++ b/deal.II/examples/step-38/doc/intro.dox @@ -0,0 +1,13 @@ + +

This will be an example of the new MeshWorker framework and will be available soon

+ +

Overview

+ +This example is devoted to the discontinuous Galerkin method, +or in short: DG method. It solves the same problem as step-12 (see +there for a description of the problem and discretization), but here +we use the MeshWorker framework in order to save reprogramming the +cell/face loops. + +

Implementation

+ diff --git a/deal.II/examples/step-38/doc/results.dox b/deal.II/examples/step-38/doc/results.dox new file mode 100644 index 0000000000..ecf8b2c67f --- /dev/null +++ b/deal.II/examples/step-38/doc/results.dox @@ -0,0 +1,91 @@ +

Results

+ + +The output of this program consist of the console output, the eps +files including the grids, and the solutions given in gnuplot format. +@code +Cycle 0: + Number of active cells: 64 + Number of degrees of freedom: 256 +Time of assemble_system1: 0.05 +Time of assemble_system2: 0.04 +solution1 and solution2 coincide. +Writing grid to ... +Writing solution to ... + +Cycle 1: + Number of active cells: 112 + Number of degrees of freedom: 448 +Time of assemble_system1: 0.09 +Time of assemble_system2: 0.07 +solution1 and solution2 coincide. +Writing grid to ... +Writing solution to ... + +Cycle 2: + Number of active cells: 214 + Number of degrees of freedom: 856 +Time of assemble_system1: 0.17 +Time of assemble_system2: 0.14 +solution1 and solution2 coincide. +Writing grid to ... +Writing solution to ... + +Cycle 3: + Number of active cells: 415 + Number of degrees of freedom: 1660 +Time of assemble_system1: 0.32 +Time of assemble_system2: 0.28 +solution1 and solution2 coincide. +Writing grid to ... +Writing solution to ... + +Cycle 4: + Number of active cells: 796 + Number of degrees of freedom: 3184 +Time of assemble_system1: 0.62 +Time of assemble_system2: 0.52 +solution1 and solution2 coincide. +Writing grid to ... +Writing solution to ... + +Cycle 5: + Number of active cells: 1561 + Number of degrees of freedom: 6244 +Time of assemble_system1: 1.23 +Time of assemble_system2: 1.03 +solution1 and solution2 coincide. +Writing grid to ... +Writing solution to ... +@endcode + +We see that, as expected, on each refinement step the two solutions +coincide. The difference measured in time of treating each face only +once (second version of the DG method) in comparison with treating +each face twice within a nested loop over all cells and all faces of +each cell (first version), is much less than one might have +expected. The gain is less than 20% on the last few refinement steps. + + + First we show the solutions on the initial mesh, the mesh after two +and after five adaptive refinement steps. + +@image html step-12.sol-0.png +@image html step-12.sol-2.png +@image html step-12.sol-5.png + + +Then we show the final grid (after 5 refinement steps) and the solution again, +this time with a nicer 3d rendering (obtained using the DataOutBase::write_vtk +function and the VTK-based VisIt visualization program) that better shows the +sharpness of the jump on the refined mesh and the over- and undershoots of the +solution along the interface: + +@image html step-12.grid-5.png +@image html step-12.3d-solution.png + + +And finally we show a plot of a 3d computation. + +@image html step-12.sol-5-3d.png + diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc new file mode 100644 index 0000000000..2e79616d83 --- /dev/null +++ b/deal.II/examples/step-38/step-38.cc @@ -0,0 +1,1650 @@ +/* $Id$ */ +/* Author: Guido Kanschat, Texas A&M University, 2009 */ + +/* $Id$ */ +/* */ +/* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 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 example 12 + // 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 + + // Here come the new include files + // for using the MeshWorker framework + + +#include +#include + + // The last step is as in all + // previous programs: +using namespace dealii; + + // @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 +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; iDGTransportEquation. Its + // member functions were already + // mentioned in the Introduction and + // will be explained + // below. Furthermore it includes + // objects of the previously defined + // Beta, RHS and + // BoundaryValues function + // classes. +template +class DGTransportEquation +{ + public: + DGTransportEquation(); + + void assemble_cell_term(const FEValues& fe_v, + FullMatrix &ui_vi_matrix, + Vector &cell_vector) const; + + void assemble_boundary_term(const FEFaceValues& fe_v, + FullMatrix &ui_vi_matrix, + Vector &cell_vector) const; + + void assemble_face_term1(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const; + + void assemble_face_term2(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_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. + // ui_vi_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 &ui_vi_matrix, + Vector &cell_vector) const +{ + // First we ask fe_v for the + // 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; pointassemble_boundary_term + // function assembles the face terms + // at boundary faces. When this + // function is invoked, fe_v is + // already reinit'ed with the current + // cell and current face. Hence it + // provides the shape values on that + // boundary face. +template +void DGTransportEquation::assemble_boundary_term( + const FEFaceValues& fe_v, + FullMatrix &ui_vi_matrix, + Vector &cell_vector) const +{ + // Again, as in the previous + // function, we ask the + // FEValues object for 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; iassemble_face_term1 + // function assembles the face terms + // corresponding to the first version + // of the DG method, cf. above. For + // that case, the face terms are + // given as a sum of integrals over + // all cell boundaries. + // + // When this function is invoked, + // fe_v and fe_v_neighbor are + // already reinit'ed with the current + // cell and the neighoring cell, + // respectively, as well as with the + // current face. Hence they provide + // the inner and outer shape values + // on the face. + // + // In addition to the cell matrix + // ui_vi_matrix this function + // gets a new argument + // ue_vi_matrix, that stores + // contributions to the system matrix + // that are based on exterior values + // of $u$ and interior values of + // $v$. Here we note that ue is + // the short notation for u + // exterior and represents $u_h^-$, + // see the introduction. +template +void DGTransportEquation::assemble_face_term1( + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const +{ + // Again, as in the previous + // function, we ask the FEValues + // objects for 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; iassemble_face_term2 function + // that assembles the face terms + // corresponding to the second + // version of the DG method, + // cf. above. For that case the face + // terms are given as a sum of + // integrals over all faces. Here we + // need two additional cell matrices + // ui_ve_matrix and + // ue_ve_matrix that will store + // contributions due to terms + // involving ui and ve as well as ue + // and ve. +template +void DGTransportEquation::assemble_face_term2( + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_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 term we've already + // seen. + for (unsigned int i=0; iDGMethod is basically + // the main class of step-6. One of + // the differences is that there's no + // ConstraintMatrix object. This is, + // because there are no hanging node + // constraints in DG discretizations. +template +class DGMethod +{ + public: + DGMethod (); + ~DGMethod (); + + void run (); + + private: + void setup_system (); + void assemble_system1 (); + void assemble_system2 (); + void solve (Vector &solution); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + const MappingQ1 mapping; + + // Furthermore we want to use DG + // elements of degree 1 (but this + // is only specified in the + // constructor). If you want to + // use a DG method of a different + // degree the whole program stays + // the same, only replace 1 in + // the constructor by the wanted + // degree. + FE_DGQ fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + // We define the quadrature + // formulae for the cell and the + // face terms of the + // discretization. + const QGauss quadrature; + const QGauss 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; + Vector solution1; + Vector solution2; + Vector right_hand_side; + + // Finally this class includes an + // object of the + // DGTransportEquations class + // described above. + const DGTransportEquation dg; +}; + + +template +DGMethod::DGMethod () + : + mapping (), + // Change here for DG + // methods of + // different degrees. + fe (1), + dof_handler (triangulation), + quadrature (4), + face_quadrature (4), + dg () +{} + + +template +DGMethod::~DGMethod () +{ + dof_handler.clear (); +} + + +template +void DGMethod::setup_system () +{ + // First we need to distribute the + // DoFs. + dof_handler.distribute_dofs (fe); + + // 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. + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + (GeometryInfo::faces_per_cell + *GeometryInfo::max_children_per_face+1)*fe.dofs_per_cell); + + // For DG discretizations we call + // the function analogue to + // DoFTools::make_sparsity_pattern. + DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern); + + // All following function calls are + // already known. + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_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 () +{ + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + std::vector dofs_neighbor (dofs_per_cell); + + // First we create the + // update_flags for the + // FEValues and the + // FEFaceValues objects. + const UpdateFlags update_flags = update_values + | update_gradients + | update_quadrature_points + | update_JxW_values; + + // Note, that on faces we do not + // need gradients but we need + // normal vectors. + const UpdateFlags face_update_flags = update_values + | update_quadrature_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. + const 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. + FEValues fe_v ( + mapping, fe, 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. + FEFaceValues fe_v_face ( + mapping, fe, face_quadrature, face_update_flags); + FESubfaceValues fe_v_subface ( + mapping, fe, face_quadrature, face_update_flags); + FEFaceValues fe_v_face_neighbor ( + mapping, fe, face_quadrature, neighbor_face_update_flags); + FESubfaceValues fe_v_subface_neighbor ( + mapping, fe, 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 vi (internal shape + // functions, i.e. shape functions + // of the current cell). To be more + // precise, the first matrix will + // include the `ui and vi terms' + // and the second will include the + // `ue and vi terms'. Here we + // recall the convention that `ui' + // is the shortcut for $u_h^+$ and + // `ue' represents $u_h^-$, see the + // introduction. + FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ue_vi_matrix (dofs_per_cell, dofs_per_cell); + + Vector cell_vector (dofs_per_cell); + + // Furthermore we need some cell + // iterators. + typename DoFHandler::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) + { + // 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 + // ui_vi_matrix and + // cell_vector to zero, + // before assembling the cell terms. + ui_vi_matrix = 0; + cell_vector = 0; + + // Now we reinit the FEValues + // object for the current cell + fe_v.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, + ui_vi_matrix, + cell_vector); + + // As in previous examples the + // vector `dofs' includes the + // dof_indices. + 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 DoFHandler::face_iterator face=cell->face(face_no); + + // and clear the + // ue_vi_matrix on each + // face. + ue_vi_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.reinit (cell, face_no); + + // and assemble the + // corresponding face + // terms. + dg.assemble_boundary_term(fe_v_face, + ui_vi_matrix, + cell_vector); + } + else + { + // Now we are not on + // the boundary of the + // domain, therefore + // there must exist a + // neighboring cell. + typename DoFHandler::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_non_children(); ++subface_no) + { + // and set the + // cell + // iterator + // neighbor_child + // to the cell + // placed + // `behind' the + // current + // subface. + typename DoFHandler::active_cell_iterator + neighbor_child + = cell->neighbor_child_on_subface (face_no, subface_no); + + Assert (!neighbor_child->has_children(), ExcInternalError()); + + // We need to + // reset the + // ue_vi_matrix + // on each + // subface + // because on + // each subface + // the un + // belong to + // different + // neighboring + // cells. + ue_vi_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.reinit (cell, face_no, subface_no); + fe_v_face_neighbor.reinit (neighbor_child, neighbor2); + + dg.assemble_face_term1(fe_v_subface, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix); + + // Then we get + // the dof + // indices of + // the + // neighbor_child + // cell + neighbor_child->get_dof_indices (dofs_neighbor); + + // and + // distribute + // ue_vi_matrix + // to the + // system_matrix + for (unsigned int i=0; iif + // (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.reinit (cell, face_no); + fe_v_face_neighbor.reinit (neighbor, neighbor2); + + dg.assemble_face_term1(fe_v_face, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_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_child_on_subface (neighbor_face_no, + neighbor_subface_no) + == cell, + ExcInternalError()); + + // Reinit the + // appropriate + // FEFaceValues + // and assemble + // the face + // terms. + fe_v_face.reinit (cell, face_no); + fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no, + neighbor_subface_no); + + dg.assemble_face_term1(fe_v_face, + fe_v_subface_neighbor, + ui_vi_matrix, + ue_vi_matrix); + } + + // Now we get the + // dof indices of + // the + // neighbor_child + // cell, + neighbor->get_dof_indices (dofs_neighbor); + + // and distribute the + // ue_vi_matrix. + for (unsigned int i=0; iface not at boundary: + } + // End of loop over all faces: + } + + // Finally we distribute the + // ui_vi_matrix + for (unsigned int i=0; iassemble_system2 function that + // implements the DG discretization + // in its second version. This + // function is very similar to the + // assemble_system1 + // function. Therefore, here we only + // discuss the differences between + // the two functions. This function + // repeatedly calls the + // assemble_face_term2 function + // of the DGTransportEquation object, + // that assembles the face terms + // written as a sum of integrals over + // all faces. Therefore, we need to + // make sure that each face is + // treated only once. This is achieved + // by introducing the rule: + // + // a) If the current and the + // neighboring cells are of the same + // refinement level we access and + // treat the face from the cell with + // lower index. + // + // b) If the two cells are of + // different refinement levels we + // access and treat the face from the + // coarser cell. + // + // Due to rule b) we do not need to + // consider case 4 (neighboring cell + // is coarser) any more. + +template +void DGMethod::assemble_system2 () +{ + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + std::vector dofs_neighbor (dofs_per_cell); + + const UpdateFlags update_flags = update_values + | update_gradients + | update_quadrature_points + | update_JxW_values; + + const UpdateFlags face_update_flags = update_values + | update_quadrature_points + | update_JxW_values + | update_normal_vectors; + + const UpdateFlags neighbor_face_update_flags = update_values; + + // Here we do not need + // fe_v_face_neighbor as case 4 + // does not occur. + FEValues fe_v ( + mapping, fe, quadrature, update_flags); + FEFaceValues fe_v_face ( + mapping, fe, face_quadrature, face_update_flags); + FESubfaceValues fe_v_subface ( + mapping, fe, face_quadrature, face_update_flags); + FEFaceValues fe_v_face_neighbor ( + mapping, fe, face_quadrature, neighbor_face_update_flags); + + + FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ue_vi_matrix (dofs_per_cell, dofs_per_cell); + + // Additionally we need the + // following two cell matrices, + // both for face term that include + // test function ve (external + // shape functions, i.e. 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 ui_ve_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ue_ve_matrix (dofs_per_cell, dofs_per_cell); + + Vector cell_vector (dofs_per_cell); + + // The following lines are roughly + // the same as in the previous + // function. + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (;cell!=endc; ++cell) + { + ui_vi_matrix = 0; + cell_vector = 0; + + fe_v.reinit (cell); + + dg.assemble_cell_term(fe_v, + ui_vi_matrix, + cell_vector); + + cell->get_dof_indices (dofs); + + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + typename DoFHandler::face_iterator face= + cell->face(face_no); + + // Case 1: + if (face->at_boundary()) + { + fe_v_face.reinit (cell, face_no); + + dg.assemble_boundary_term(fe_v_face, + ui_vi_matrix, + cell_vector); + } + else + { + Assert (cell->neighbor(face_no).state() == IteratorState::valid, + ExcInternalError()); + typename DoFHandler::cell_iterator neighbor= + cell->neighbor(face_no); + // Case 2: + if (face->has_children()) + { + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face_no); + + for (unsigned int subface_no=0; + subface_non_children(); ++subface_no) + { + typename DoFHandler::cell_iterator neighbor_child + = cell->neighbor_child_on_subface (face_no, subface_no); + Assert (!neighbor_child->has_children(), ExcInternalError()); + + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; + + fe_v_subface.reinit (cell, face_no, subface_no); + fe_v_face_neighbor.reinit (neighbor_child, neighbor2); + + dg.assemble_face_term2(fe_v_subface, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); + + 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); + + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; + + fe_v_face.reinit (cell, face_no); + fe_v_face_neighbor.reinit (neighbor, neighbor2); + + dg.assemble_face_term2(fe_v_face, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); + + neighbor->get_dof_indices (dofs_neighbor); + + for (unsigned int i=0; i +void DGMethod::solve (Vector &solution) +{ + SolverControl solver_control (1000, 1e-12, false, false); + SolverRichardson<> solver (solver_control); + + // Here we create the + // preconditioner, + PreconditionBlockSSOR > preconditioner; + + // we assigned the matrix to it and + // set the right block size. + preconditioner.initialize(system_matrix, fe.dofs_per_cell); + + // After these preparations we are + // ready to start the linear solver. + solver.solve (system_matrix, solution, right_hand_side, + preconditioner); +} + + + // 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 DoFHandler::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); + + 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"; + 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 (); + + data_out.write_gnuplot(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<6; ++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 () +{ + try + { + DGMethod<2> dgmethod; + dgmethod.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; +} + +