From: heister Date: Thu, 3 Jan 2013 21:37:23 +0000 (+0000) Subject: reintegrate branch merge_mg_into_dof_handler X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f8a1399a13bffc57d413c701a2a412423686c7dd;p=dealii-svn.git reintegrate branch merge_mg_into_dof_handler git-svn-id: https://svn.dealii.org/trunk@27912 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/common/dealiitemplates.pm.in b/deal.II/common/dealiitemplates.pm.in index 80b7a15838..4642d80817 100644 --- a/deal.II/common/dealiitemplates.pm.in +++ b/deal.II/common/dealiitemplates.pm.in @@ -102,7 +102,7 @@ our @dimensions = (1,2,3); our @codim1_dimensions = (1,2); # DoFHandlers -our @dofhandler_templates = qw(DoFHandler hp::DoFHandler MGDoFHandler); +our @dofhandler_templates = qw(DoFHandler hp::DoFHandler); our @dofhandlers = map { $t = $_ ; map { "$_<$t>" } @dofhandler_templates } @dimensions; # All possible number types diff --git a/deal.II/common/template-arguments.in b/deal.II/common/template-arguments.in index 09d4d3354d..8178ad9af6 100644 --- a/deal.II/common/template-arguments.in +++ b/deal.II/common/template-arguments.in @@ -1,3 +1,4 @@ +BOOL := { true; false } REAL_SCALARS := { double; float; long double } COMPLEX_SCALARS := { std::complex; std::complex; @@ -36,23 +37,31 @@ SERIAL_VECTORS := { Vector; @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR@; } +VECTORS_WITH_MATRIX := { Vector; + Vector ; + Vector; + + BlockVector; + BlockVector; + BlockVector; + + @DEAL_II_EXPAND_TRILINOS_VECTOR@; + @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; + } + DOFHANDLERS := { DoFHandler; - hp::DoFHandler; - MGDoFHandler } + hp::DoFHandler } DOFHANDLER_TEMPLATES := { DoFHandler; - hp::DoFHandler; - MGDoFHandler } + hp::DoFHandler } TRIANGULATION_AND_DOFHANDLER_TEMPLATES := { Triangulation; DoFHandler; - hp::DoFHandler; - MGDoFHandler} + hp::DoFHandler } TRIANGULATION_AND_DOFHANDLERS := { Triangulation; DoFHandler; - hp::DoFHandler; - MGDoFHandler } + hp::DoFHandler } FEVALUES_BASES := { FEValuesBase; diff --git a/deal.II/doc/doxygen/headers/iterators.h b/deal.II/doc/doxygen/headers/iterators.h index 0ab81cc360..52f44bea81 100644 --- a/deal.II/doc/doxygen/headers/iterators.h +++ b/deal.II/doc/doxygen/headers/iterators.h @@ -239,22 +239,19 @@ of accessor classes: topological connection of a mesh is available from an accessor such as to request iterators pointing to neighbors of a cell. -- The DoFAccessor class lets you access information related to degrees of - freedom associated with cells, faces, etc; it does so for both DoFHandler - and hp::DoFHandler objects. Note that the DoFAccessor class is derived from - either TriaAccessor or CellAccessor (depending on whether the DoFAccessor - points to an object of full dimension or not) and so is able to provide a - superset of information over its base classes. +- The DoFAccessor class lets you access information related to degrees + of freedom associated with cells, faces, etc; it does so for both + DoFHandler and hp::DoFHandler objects. Note that the DoFAccessor + class is derived from either TriaAccessor or CellAccessor (depending + on whether the DoFAccessor points to an object of full dimension or + not) and so is able to provide a superset of information over its + base classes. Additionally, the DoFAccessor class comes in two + flavors, one accessing degrees of freedom on the level of a cell and + the other accessing the active dofs of an active cell. - The DoFCellAccessor class has the same purpose and relation to DoFCellAccessor as the CellAccessor has to TriaAccessor. -- The MGDoFAccessor and MFCellDoFAccessor give access to the multilevel degree - of freedom information associated with MGDoFHandler objects. As before, - since the MGDoFHandler extends the data structures of regular DoFHandler - objects, the MGDoFAccessor is also derived from the DoFAccessor (or - DoFCellAccessor) class to provide a superset of information. - Except to look up member documentation, you will not usually have to deal with the actual class names listed above. Rather, one uses the typedefs provided by the container classes Triangulation, DoFHandler, hp::DoFHandler and diff --git a/deal.II/examples/step-50/Makefile b/deal.II/examples/step-50/Makefile new file mode 100644 index 0000000000..3af5d515e1 --- /dev/null +++ b/deal.II/examples/step-50/Makefile @@ -0,0 +1,143 @@ +# $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 = step-50 + +# 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 *vtk *ucd *.d2 + + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# The next statement tells 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. deal.II has two +# libraries: one for the debug mode version of the +# application and one for optimized mode. +libs.g := $(lib-deal2.g) +libs.o := $(lib-deal2.o) + + +# We now use the variable defined above to 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 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. +all: $(target)$(EXEEXT) +$(target)$(EXEEXT) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@ $^ $(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)$(EXEEXT) + @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========= $( $@" + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +./%.$(OBJEXT) : + @echo "==============optimized===== $( $@" + @$(CXX) $(CXXFLAGS.o) -c $< -o $@ + + +# The following statement tells make that the rules `run' and `clean' +# are not expected to produce files of the same name as Makefile rules +# usually do. +.PHONY: all run clean + + +# Finally there is a rule which you normally need not care much about: +# since the executable depends on some include files from the library, +# besides the C++ application file of course, it is necessary to +# re-generate the executable when one of the files it depends on has +# changed. The following rule creates a dependency file +# `Makefile.dep', which `make' uses to determine when to regenerate +# the executable. This file is automagically remade whenever needed, +# i.e. whenever one of the cc-/h-files changed. Make detects whether +# to remake this file upon inclusion at the bottom of this file. +# +# If the creation of Makefile.dep fails, blow it away and fail +Makefile.dep: $(target).cc Makefile \ + $(shell echo $D/include/deal.II/*/*.h) + @echo ============================ Remaking $@ + @$D/common/scripts/make_dependencies $(INCLUDE) -B. $(target).cc \ + > $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then true ; else rm $@ ; false ; fi + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-50/doc/builds-on b/deal.II/examples/step-50/doc/builds-on new file mode 100644 index 0000000000..17402734c7 --- /dev/null +++ b/deal.II/examples/step-50/doc/builds-on @@ -0,0 +1 @@ +step-6 diff --git a/deal.II/examples/step-50/doc/intro.dox b/deal.II/examples/step-50/doc/intro.dox new file mode 100644 index 0000000000..99666726d3 --- /dev/null +++ b/deal.II/examples/step-50/doc/intro.dox @@ -0,0 +1,90 @@ +
+ +This program has evolved from a version originally written by Guido +Kanschat in 2003. It has undergone significant revisions by Bärbel +Janssen, Guido Kanschat and Wolfgang Bangerth in 2009 and 2010 to demonstrate +multigrid algorithms on adaptively refined meshes. + + + + +

Introduction

+ + +This example shows the basic usage of the multilevel functions in +deal.II. It solves the same problem as used in step-6, +but demonstrating the things one has to provide when using multigrid +as a preconditioner. In particular, this requires that we define a +hierarchy of levels, provide transfer operators from one level to the +next and back, and provide representations of the Laplace operator on +each level. + +In order to allow sufficient flexibility in conjunction with systems of +differential equations and block preconditioners, quite a few different objects +have to be created before starting the multilevel method, although +most of what needs to be done is provided by deal.II itself. These are +
    +
  • An the object handling transfer between grids; we use the + MGTransferPrebuilt class for this that does almost all of the work + inside the library. +
  • The solver on the coarsest level; here, we use MGCoarseGridHouseholder. +
  • The smoother on all other levels, which in our case will be the + MGSmootherRelaxation class using SOR as the underlying method +
  • And MGMatrix, a class having a special level multiplication, i.e. we + basically store one matrix per grid level and allow multiplication + with it. +
+Most of these objects will only be needed inside the function that +actually solves the linear system. There, these objects are combined +in an object of type Multigrid, containing the implementation of the +V-cycle, which is in turn used by the preconditioner PreconditionMG, +ready for plug-in into a linear solver of the LAC library. + +The multilevel method in deal.II follows in many respects the outlines +of the various publications by James Bramble, Joseph Pasciak and +Jinchao Xu (i.e. the "BPX" framework). In order to understand many of +the options, a rough familiarity with their work is quite helpful. + +However, in comparison to this framework, the implementation in +deal.II has to take into account the fact that we want to solve linear +systems on adaptively refined meshes. This leads to the complication +that it isn't quite as clear any more what exactly a "level" in a +multilevel hierarchy of a mesh is. The following image shows what we +consider to be a "level": + +

+ @image html "hanging_nodes.png" "" +

+ +In other words, the fine level in this mesh consists only of the +degrees of freedom that are defined on the refined cells, but does not +extend to that part of the domain that is not refined. While this +guarantees that the overall effort grows as ${\cal O}(N)$ as necessary +for optimal multigrid complexity, it leads to problems when defining +where to smooth and what boundary conditions to pose for the operators +defined on individual levels if the level boundary is not an external +boundary. These questions are discussed in detail in the +@ref mg_paper "Multigrid paper by Janssen and Kanschat" that describes +the implementation in deal.II. + + + +

The testcase

+ +The problem we solve here is exactly the same as in +step-6, the only difference being the solver we use +here. You may want to look there for a definition of what we solve, +right hand side and boundary conditions. Obviously, the program would +also work if we changed the geometry and other pieces of data that +defines this particular problem. + +The things that are new are all those parts that concern the +multigrid. In particular, this includes the following members of the +main class: +- LaplaceProblem::mg_dof_handler +- LaplaceProblem::mg_sparsity +- LaplaceProblem::mg_matrices +- LaplaceProblem::mg_interface_matrices_up +- LaplaceProblem::assemble_multigrid () +- LaplaceProblem::solve () +Take a look at these functions. diff --git a/deal.II/examples/step-50/doc/kind b/deal.II/examples/step-50/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/deal.II/examples/step-50/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/deal.II/examples/step-50/doc/results.dox b/deal.II/examples/step-50/doc/results.dox new file mode 100644 index 0000000000..762f8091b4 --- /dev/null +++ b/deal.II/examples/step-50/doc/results.dox @@ -0,0 +1,101 @@ +

Results

+ +The output that this program generates is, of course, the same as that +of step-6, so you may see there for more results. On the +other hand, since no tutorial program is a good one unless it has at +least one colorful picture, here is, again, the solution: + +

+ @image html step-16.solution.png +

+ +When run, the output of this program is +
+Cycle 0:
+   Number of active cells:       20
+   Number of degrees of freedom: 25 (by level: 8, 25)
+   7 CG iterations needed to obtain convergence.
+Cycle 1:
+   Number of active cells:       44
+   Number of degrees of freedom: 57 (by level: 8, 25, 48)
+   8 CG iterations needed to obtain convergence.
+Cycle 2:
+   Number of active cells:       92
+   Number of degrees of freedom: 117 (by level: 8, 25, 80, 60)
+   9 CG iterations needed to obtain convergence.
+Cycle 3:
+   Number of active cells:       188
+   Number of degrees of freedom: 221 (by level: 8, 25, 80, 200)
+   12 CG iterations needed to obtain convergence.
+Cycle 4:
+   Number of active cells:       416
+   Number of degrees of freedom: 485 (by level: 8, 25, 89, 288, 280)
+   13 CG iterations needed to obtain convergence.
+Cycle 5:
+   Number of active cells:       800
+   Number of degrees of freedom: 925 (by level: 8, 25, 89, 288, 784, 132)
+   14 CG iterations needed to obtain convergence.
+Cycle 6:
+   Number of active cells:       1628
+   Number of degrees of freedom: 1865 (by level: 8, 25, 89, 304, 1000, 1164, 72)
+   14 CG iterations needed to obtain convergence.
+Cycle 7:
+   Number of active cells:       3194
+   Number of degrees of freedom: 3603 (by level: 8, 25, 89, 328, 1032, 2200, 1392)
+   16 CG iterations needed to obtain convergence.
+
+That's not perfect — we would have hoped for a constant number +of iterations rather than one that increases as we get more and more +degrees of freedom — but it is also not far away. The reason for +this is easy enough to understand, however: since we have a strongly +varying coefficient, the operators that we assembly by quadrature on +the lower levels become worse and worse approximations of the operator +on the finest level. Consequently, even if we had perfect solvers on +the coarser levels, they would not be good preconditioners on the +finest level. This theory is easily tested by comparing results when +we use a constant coefficient: in that case, the number of iterations +remains constant at 9 after the first three or four refinement steps. + +We can also compare what this program produces with how @ref step_5 +"step-5" performed. To solve the same problem as in step-5, the only +two changes that are necessary are (i) to replace the body of the +function LaplaceProblem::refine_grid by a call to +triangulation.refine_global(1), and (ii) to use the same +SolverControl object and tolerance as in step-5 — the rest of the +program remains unchanged. In that case, here is how the solvers used +in step-5 and the multigrid solver used in the current program +compare: + + + + + + + + +
cellsstep-5step-16
20 13 6
80 17 7
320 29 9
1280 51 10
5120 94 11
2048018013
+This isn't only fewer iterations than in step-5 (each of which +is, however, much more expensive) but more importantly, the number of +iterations also grows much more slowly under mesh refinement (again, +it would be almost constant if the coefficient was constant rather +than strongly varying as chosen here). This justifies the common +observation that, whenever possible, multigrid methods should be used +for second order problems. + + +

Possible extensions

+ +A close inspection of this program's performance shows that it is mostly +dominated by matrix-vector operations. step-37 shows one way +how this can be avoided by working with matrix-free methods. + +Another avenue would be to use algebraic multigrid methods. The +geometric multigrid method used here can at times be a bit awkward to +implement because it needs all those additional data structures, and +it becomes even more difficult if the program is to run in %parallel on +machines coupled through MPI, for example. In that case, it would be +simpler if one could use a black-box preconditioner that uses some +sort of multigrid hierarchy for good performance but can figure out +level matrices and similar things out by itself. Algebraic multigrid +methods do exactly this, and we will use them in +step-31 for the solution of a Stokes problem. diff --git a/deal.II/examples/step-50/doc/step-16.solution.png b/deal.II/examples/step-50/doc/step-16.solution.png new file mode 100644 index 0000000000..2550b85a26 Binary files /dev/null and b/deal.II/examples/step-50/doc/step-16.solution.png differ diff --git a/deal.II/examples/step-50/doc/tooltip b/deal.II/examples/step-50/doc/tooltip new file mode 100644 index 0000000000..2fd65590a0 --- /dev/null +++ b/deal.II/examples/step-50/doc/tooltip @@ -0,0 +1 @@ +Multigrid on adaptive meshes. diff --git a/deal.II/examples/step-50/step-50.cc b/deal.II/examples/step-50/step-50.cc new file mode 100644 index 0000000000..8cc10afd97 --- /dev/null +++ b/deal.II/examples/step-50/step-50.cc @@ -0,0 +1,1064 @@ +/* Author: Guido Kanschat, University of Heidelberg, 2003 */ +/* Baerbel Janssen, University of Heidelberg, 2010 */ +/* Wolfgang Bangerth, Texas A&M University, 2010 */ + +/* $Id$ */ +/* */ +/* Copyright (C) 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012 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. */ + +// As discussed in the introduction, most of +// this program is copied almost verbatim +// from step-6, which itself is only a slight +// modification of step-5. Consequently, a +// significant part of this program is not +// new if you've read all the material up to +// step-6, and we won't comment on that part +// of the functionality that is +// unchanged. Rather, we will focus on those +// aspects of the program that have to do +// with the multigrid functionality which +// forms the new aspect of this tutorial +// program. + +// @sect3{Include files} + +// Again, the first few include files +// are already known, so we won't +// comment on them: +#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 + +#include +#include +#include + +// These, now, are the include necessary for +// the multi-level methods. The first two +// declare classes that allow us to enumerate +// degrees of freedom not only on the finest +// mesh level, but also on intermediate +// levels (that's what the MGDoFHandler class +// does) as well as allow to access this +// information (iterators and accessors over +// these cells). +// +// The rest of the include files deals with +// the mechanics of multigrid as a linear +// operator (solver or preconditioner). +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// This is C++: +#include +#include + +// The last step is as in all +// previous programs: +namespace Step16 +{ + using namespace dealii; + + + // @sect3{The LaplaceProblem class template} + + // This main class is basically the same + // class as in step-6. As far as member + // functions is concerned, the only addition + // is the assemble_multigrid + // function that assembles the matrices that + // correspond to the discrete operators on + // intermediate levels: + template + class LaplaceProblem + { + public: + LaplaceProblem (const unsigned int deg); + void run (); + + private: + void setup_system (); + void assemble_system (); + void assemble_multigrid (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + parallel::distributed::Triangulation triangulation; + FE_Q fe; + DoFHandler mg_dof_handler; + + typedef TrilinosWrappers::SparseMatrix matrix_t; + typedef TrilinosWrappers::MPI::Vector vector_t; + + SparsityPattern sparsity_pattern; + matrix_t system_matrix; + + // We need an additional object for the + // hanging nodes constraints. They are + // handed to the transfer object in the + // multigrid. Since we call a compress + // inside the multigrid these constraints + // are not allowed to be inhomogeneous so + // we store them in different ConstraintMatrix + // objects. + ConstraintMatrix hanging_node_constraints; + ConstraintMatrix constraints; + + vector_t solution; + vector_t system_rhs; + + const unsigned int degree; + + // The following four objects are the + // only additional member variables, + // compared to step-6. They first three + // represent the + // operators that act on individual + // levels of the multilevel hierarchy, + // rather than on the finest mesh as do + // the objects above while the last object + // stores information about the boundary + // indices on each level and information + // about indices lying on a refinement + // edge between two different refinement + // levels. + // + // To facilitate having objects on each + // level of a multilevel hierarchy, + // deal.II has the MGLevelObject class + // template that provides storage for + // objects on each level. What we need + // here are matrices on each level, which + // implies that we also need sparsity + // patterns on each level. As outlined in + // the @ref mg_paper, the operators + // (matrices) that we need are actually + // twofold: one on the interior of each + // level, and one at the interface + // between each level and that part of + // the domain where the mesh is + // coarser. In fact, we will need the + // latter in two versions: for the + // direction from coarse to fine mesh and + // from fine to coarse. Fortunately, + // however, we here have a self-adjoint + // problem for which one of these is the + // transpose of the other, and so we only + // have to build one; we choose the one + // from coarse to fine. + MGLevelObject mg_sparsity_patterns; + MGLevelObject mg_matrices; + MGLevelObject mg_interface_matrices; + MGConstrainedDoFs mg_constrained_dofs; + }; + + + + // @sect3{Nonconstant coefficients} + + // The implementation of nonconstant + // coefficients is copied verbatim + // from step-5 and step-6: + + 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; i::value (points[i]); + } + + + // @sect3{The LaplaceProblem class implementation} + + // @sect4{LaplaceProblem::LaplaceProblem} + + // The constructor is left mostly + // unchanged. We take the polynomial degree + // of the finite elements to be used as a + // constructor argument and store it in a + // member variable. + // + // By convention, all adaptively refined + // triangulations in deal.II never change by + // more than one level across a face between + // cells. For our multigrid algorithms, + // however, we need a slightly stricter + // guarantee, namely that the mesh also does + // not change by more than refinement level + // across vertices that might connect two + // cells. In other words, we must prevent the + // following situation: + // + // @image html limit_level_difference_at_vertices.png "" + // + // This is achieved by passing the + // Triangulation::limit_level_difference_at_vertices + // flag to the constructor of the + // triangulation class. + template + LaplaceProblem::LaplaceProblem (const unsigned int degree) + : + triangulation (MPI_COMM_WORLD,Triangulation:: + limit_level_difference_at_vertices), + fe (degree), + mg_dof_handler (triangulation), + degree(degree) + {} + + + + // @sect4{LaplaceProblem::setup_system} + + // The following function extends what the + // corresponding one in step-6 did. The top + // part, apart from the additional output, + // does the same: + template + void LaplaceProblem::setup_system () + { + mg_dof_handler.distribute_dofs (fe); + mg_dof_handler.distribute_mg_dofs (fe); + + + // Here we output not only the + // degrees of freedom on the finest + // level, but also in the + // multilevel structure + deallog << "Number of degrees of freedom: " + << mg_dof_handler.n_dofs(); + + for (unsigned int l=0; l::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (1); + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + VectorTools::interpolate_boundary_values (static_cast&>(mg_dof_handler), + dirichlet_boundary, + constraints); + constraints.close (); + hanging_node_constraints.close (); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + system_matrix.reinit (sparsity_pattern); + + // The multigrid constraints have to be + // initialized. They need to know about + // the boundary values as well, so we + // pass the dirichlet_boundary + // here as well. + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); + + + // Now for the things that concern the + // multigrid data structures. First, we + // resize the multi-level objects to hold + // matrices and sparsity patterns for every + // level. The coarse level is zero (this is + // mandatory right now but may change in a + // future revision). Note that these + // functions take a complete, inclusive + // range here (not a starting index and + // size), so the finest level is + // n_levels-1. We first have + // to resize the container holding the + // SparseMatrix classes, since they have to + // release their SparsityPattern before the + // can be destroyed upon resizing. + const unsigned int n_levels = triangulation.n_levels(); + + mg_interface_matrices.resize(0, n_levels-1); + mg_interface_matrices.clear (); + mg_matrices.resize(0, n_levels-1); + mg_matrices.clear (); + mg_sparsity_patterns.resize(0, n_levels-1); + + // Now, we have to provide a matrix on each + // level. To this end, we first use the + // MGTools::make_sparsity_pattern function + // to first generate a preliminary + // compressed sparsity pattern on each + // level (see the @ref Sparsity module for + // more information on this topic) and then + // copy it over to the one we really + // want. The next step is to initialize + // both kinds of level matrices with these + // sparsity patterns. + // + // It may be worth pointing out that the + // interface matrices only have entries for + // degrees of freedom that sit at or next + // to the interface between coarser and + // finer levels of the mesh. They are + // therefore even sparser than the matrices + // on the individual levels of our + // multigrid hierarchy. If we were more + // concerned about memory usage (and + // possibly the speed with which we can + // multiply with these matrices), we should + // use separate and different sparsity + // patterns for these two kinds of + // matrices. + for (unsigned int level=0; level + void LaplaceProblem::assemble_system () + { + const QGauss quadrature_formula(degree+1); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + 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); + } + } + + + // @sect4{LaplaceProblem::assemble_multigrid} + + // The next function is the one that builds + // the linear operators (matrices) that + // define the multigrid method on each level + // of the mesh. The integration core is the + // same as above, but the loop below will go + // over all existing cells instead of just + // the active ones, and the results must be + // entered into the correct matrix. Note also + // that since we only do multi-level + // preconditioning, no right-hand side needs + // to be assembled here. + // + // Before we go there, however, we have to + // take care of a significant amount of book + // keeping: + template + void LaplaceProblem::assemble_multigrid () + { + QGauss quadrature_formula(1+degree); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + // Next a few things that are specific to + // building the multigrid data structures + // (since we only need them in the current + // function, rather than also elsewhere, we + // build them here instead of the + // setup_system + // function). Some of the following may be + // a bit obscure if you're not familiar + // with the algorithm actually implemented + // in deal.II to support multilevel + // algorithms on adaptive meshes; if some + // of the things below seem strange, take a + // look at the @ref mg_paper. + // + // Our first job is to identify those + // degrees of freedom on each level that + // are located on interfaces between + // adaptively refined levels, and those + // that lie on the interface but also on + // the exterior boundary of the domain. As + // in many other parts of the library, we + // do this by using boolean masks, + // i.e. vectors of booleans each element of + // which indicates whether the + // corresponding degree of freedom index is + // an interface DoF or not. The MGConstraints + // already computed the information for us + // when we called initialize in setup_system(). + std::vector > interface_dofs + = mg_constrained_dofs.get_refinement_edge_indices (); + std::vector > boundary_interface_dofs + = mg_constrained_dofs.get_refinement_edge_boundary_indices (); + + // The indices just identified will later + // be used to decide where the assembled value + // has to be added into on each level. + // On the other hand, + // we also have to impose zero boundary + // conditions on the external boundary of + // each level. But this the MGConstraints + // knows it. So we simply ask for them by calling + // get_boundary_indices (). + // The third step is to construct + // constraints on all those degrees of + // freedom: their value should be zero + // after each application of the level + // operators. To this end, we construct + // ConstraintMatrix objects for each level, + // and add to each of these constraints for + // each degree of freedom. Due to the way + // the ConstraintMatrix stores its data, + // the function to add a constraint on a + // single degree of freedom and force it to + // be zero is called + // Constraintmatrix::add_line(); doing so + // for several degrees of freedom at once + // can be done using + // Constraintmatrix::add_lines(): + std::vector boundary_constraints (triangulation.n_levels()); + std::vector boundary_interface_constraints (triangulation.n_levels()); + for (unsigned int level=0; levelassemble_system, with two + // exceptions: (i) we don't need a right + // hand side, and more significantly (ii) we + // don't just loop over all active cells, + // but in fact all cells, active or + // not. Consequently, the correct iterator + // to use is MGDoFHandler::cell_iterator + // rather than + // MGDoFHandler::active_cell_iterator. Let's + // go about it: + typename DoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); + + for (; cell!=endc; ++cell) + if (cell->level_subdomain_id()==triangulation.locally_owned_subdomain()) + { + cell_matrix = 0; + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + + // Next, we need to copy local + // contributions into the level + // objects. We can do this in the same + // way as in the global assembly, using + // a constraint object that takes care + // of constrained degrees (which here + // are only boundary nodes, as the + // individual levels have no hanging + // node constraints). Note that the + // boundary_constraints + // object makes sure that the level + // matrices contains no contributions + // from degrees of freedom at the + // interface between cells of different + // refinement level. + boundary_constraints[cell->level()] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_matrices[cell->level()]); + + // The next step is again slightly more + // obscure (but explained in the @ref + // mg_paper): We need the remainder of + // the operator that we just copied + // into the mg_matrices + // object, namely the part on the + // interface between cells at the + // current level and cells one level + // coarser. This matrix exists in two + // directions: for interior DoFs (index + // $i$) of the current level to those + // sitting on the interface (index + // $j$), and the other way around. Of + // course, since we have a symmetric + // operator, one of these matrices is + // the transpose of the other. + // + // The way we assemble these matrices + // is as follows: since the are formed + // from parts of the local + // contributions, we first delete all + // those parts of the local + // contributions that we are not + // interested in, namely all those + // elements of the local matrix for + // which not $i$ is an interface DoF + // and $j$ is not. The result is one of + // the two matrices that we are + // interested in, and we then copy it + // into the + // mg_interface_matrices + // object. The + // boundary_interface_constraints + // object at the same time makes sure + // that we delete contributions from + // all degrees of freedom that are not + // only on the interface but also on + // the external boundary of the domain. + // + // The last part to remember is how to + // get the other matrix. Since it is + // only the transpose, we will later + // (in the solve() + // function) be able to just pass the + // transpose matrix where necessary. + for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && + interface_dofs[cell->level()][local_dof_indices[j]]==false)) + cell_matrix(i,j) = 0; + + boundary_interface_constraints[cell->level()] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_interface_matrices[cell->level()]); + } + } + + + + // @sect4{LaplaceProblem::solve} + + // This is the other function that is + // significantly different in support of the + // multigrid solver (or, in fact, the + // preconditioner for which we use the + // multigrid method). + // + // Let us start out by setting up two of the + // components of multilevel methods: transfer + // operators between levels, and a solver on + // the coarsest level. In finite element + // methods, the transfer operators are + // derived from the finite element function + // spaces involved and can often be computed + // in a generic way independent of the + // problem under consideration. In that case, + // we can use the MGTransferPrebuilt class + // that, given the constraints on the global + // level and an MGDoFHandler object computes + // the matrices corresponding to these + // transfer operators. + // + // The second part of the following lines + // deals with the coarse grid solver. Since + // our coarse grid is very coarse indeed, we + // decide for a direct solver (a Householder + // decomposition of the coarsest level + // matrix), even if its implementation is not + // particularly sophisticated. If our coarse + // mesh had many more cells than the five we + // have here, something better suited would + // obviously be necessary here. + template + void LaplaceProblem::solve () + { + + // Create the object that deals with the transfer + // between different refinement levels. We need to + // pass it the hanging node constraints. + MGTransferPrebuilt mg_transfer(hanging_node_constraints, mg_constrained_dofs); + // Now the prolongation matrix has to be built. + // This matrix needs to take the boundary values on + // each level into account and needs to know about + // the indices at the refinement egdes. The + // MGConstraints knows about that so + // pass it as an argument. + mg_transfer.build_matrices(mg_dof_handler); + + matrix_t coarse_matrix; + coarse_matrix.copy_from (mg_matrices[0]); + //MGCoarseGridHouseholder coarse_grid_solver; + //coarse_grid_solver.initialize (coarse_matrix); + + SolverControl coarse_solver_control (1000, 1e-10, false, false); + SolverGMRES coarse_solver(coarse_solver_control); + PreconditionIdentity id; + MGCoarseGridLACIteration,vector_t> coarse_grid_solver(coarse_solver, + coarse_matrix, + id); + + // The next component of a multilevel + // solver or preconditioner is that we need + // a smoother on each level. A common + // choice for this is to use the + // application of a relaxation method (such + // as the SOR, Jacobi or Richardson method). The + // MGSmootherPrecondition class provides + // support for this kind of + // smoother. Here, we opt for the + // application of a single SOR + // iteration. To this end, we define an + // appropriate typedef and + // then setup a smoother object. + // + // Since this smoother needs temporary + // vectors to store intermediate results, + // we need to provide a VectorMemory + // object. Since these vectors will be + // reused over and over, the + // GrowingVectorMemory is more time + // efficient than the PrimitiveVectorMemory + // class in the current case. + // + // The last step is to initialize the + // smoother object with our level matrices + // and to set some smoothing parameters. + // The initialize() function + // can optionally take additional arguments + // that will be passed to the smoother + // object on each level. In the current + // case for the SOR smoother, this could, + // for example, include a relaxation + // parameter. However, we here leave these + // at their default values. The call to + // set_steps() indicates that + // we will use two pre- and two + // post-smoothing steps on each level; to + // use a variable number of smoother steps + // on different levels, more options can be + // set in the constructor call to the + // mg_smoother object. + // + // The last step results from the fact that + // we use the SOR method as a smoother - + // which is not symmetric - but we use the + // conjugate gradient iteration (which + // requires a symmetric preconditioner) + // below, we need to let the multilevel + // preconditioner make sure that we get a + // symmetric operator even for nonsymmetric + // smoothers: + typedef TrilinosWrappers::PreconditionSOR Smoother; + GrowingVectorMemory vector_memory; + MGSmootherPrecondition + mg_smoother(vector_memory); + mg_smoother.initialize(mg_matrices); + mg_smoother.set_steps(2); + mg_smoother.set_symmetric(true); + + // The next preparatory step is that we + // must wrap our level and interface + // matrices in an object having the + // required multiplication functions. We + // will create two objects for the + // interface objects going from coarse to + // fine and the other way around; the + // multigrid algorithm will later use the + // transpose operator for the latter + // operation, allowing us to initialize + // both up and down versions of the + // operator with the matrices we already + // built: + MGMatrix mg_matrix(&mg_matrices); + MGMatrix mg_interface_up(&mg_interface_matrices); + MGMatrix mg_interface_down(&mg_interface_matrices); + + // Now, we are ready to set up the + // V-cycle operator and the + // multilevel preconditioner. + Multigrid mg(mg_dof_handler, + mg_matrix, + coarse_grid_solver, + mg_transfer, + mg_smoother, + mg_smoother); + mg.set_edge_matrices(mg_interface_down, mg_interface_up); + + PreconditionMG > + preconditioner(mg_dof_handler, mg, mg_transfer); + + // With all this together, we can finally + // get about solving the linear system in + // the usual way: + SolverControl solver_control (1000, 1e-12); + SolverCG cg (solver_control); + + solution = 0; + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + constraints.distribute (solution); + } + + + + // @sect4{Postprocessing} + + // The following two functions postprocess a + // solution once it is computed. In + // particular, the first one refines the mesh + // at the beginning of each cycle while the + // second one outputs results at the end of + // each such cycle. The functions are almost + // unchanged from those in step-6, with the + // exception of two minor differences: The + // KellyErrorEstimator::estimate function + // wants an argument of type DoFHandler, not + // MGDoFHandler, and so we have to cast from + // derived to base class; and we generate + // output in VTK format, to use the more + // modern visualization programs available + // today compared to those that were + // available when step-6 was written. + template + void LaplaceProblem::refine_grid () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (static_cast&>(mg_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.03); + triangulation.execute_coarsening_and_refinement (); + } + + + + template + void LaplaceProblem::output_results (const unsigned int cycle) const + { + DataOut data_out; + + data_out.attach_dof_handler (mg_dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + std::ostringstream filename; + filename << "solution-" + << cycle + << ".vtk"; + + std::ofstream output (filename.str().c_str()); + data_out.write_vtk (output); + } + + + // @sect4{LaplaceProblem::run} + + // Like several of the functions above, this + // is almost exactly a copy of of the + // corresponding function in step-6. The only + // difference is the call to + // assemble_multigrid that takes + // care of forming the matrices on every + // level that we need in the multigrid + // method. + template + void LaplaceProblem::run () + { + for (unsigned int cycle=0; cycle<8; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + triangulation.refine_global (1); + } + else + refine_grid (); + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << mg_dof_handler.n_dofs() + << " (by level: "; + for (unsigned int level=0; level laplace_problem(1); + laplace_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/deal.II/include/deal.II/dofs/block_info.h b/deal.II/include/deal.II/dofs/block_info.h index 9738685202..835d5270e5 100644 --- a/deal.II/include/deal.II/dofs/block_info.h +++ b/deal.II/include/deal.II/dofs/block_info.h @@ -23,7 +23,6 @@ DEAL_II_NAMESPACE_OPEN // Forward declarations template class DoFHandler; -template class MGDoFHandler; namespace hp { template class DoFHandler; @@ -48,7 +47,7 @@ namespace hp * BlockIndices::block_start() will return the start index for each of * the blocks. * - * When an MGDoFHandler is used, the same structure is automatically + * When a DoFHandler with levels is used, the same structure is automatically * generated for each level. The level blocks can be accessed through * level(). * @@ -94,26 +93,25 @@ public: * describing block structure * of the DoFHandler. * + * By default, this function will + * attempt to initialize whatever + * is possible. If active dofs + * have been assigned int the + * DoFHandler argument, they + * BlockIndices for those will be + * generated. The same for level + * dofs. + * + * This default behavior can be + * overridden by the two + * parameters, which can switch + * off active dofs or level dofs. + * * This function will also clear * the local() indices. */ template - void initialize(const DoFHandler &); - - /** - * @brief Fill the object with values - * describing level block - * structure of the - * MGDoFHandler. If - * levels_only is false, - * the other initialize() is - * called as well. - * - * This function will also clear - * the local() indices. - */ - template - void initialize(const MGDoFHandler &, bool levels_only = false); + void initialize(const DoFHandler &, bool levels_only = false, bool active_only = false); /** * @brief Initialize block structure diff --git a/deal.II/include/deal.II/dofs/dof_accessor.h b/deal.II/include/deal.II/dofs/dof_accessor.h index d587001c98..beebb35966 100644 --- a/deal.II/include/deal.II/dofs/dof_accessor.h +++ b/deal.II/include/deal.II/dofs/dof_accessor.h @@ -120,32 +120,25 @@ namespace internal /** - * A class that gives access to the degrees of freedom stored in a DoFHandler - * or hp::DoFHandler object. Accessors are used to, well, access the data that + * A class that gives access to the degrees of freedom stored in a DoFHandler, + * MGSoDHandler, or hp::DoFHandler object. Accessors are used to, access the data that * pertains to edges, faces, and cells of a triangulation. The concept is * explained in more detail in connection to @ref Iterators. * - * This class follows mainly the route laid out by the accessor library - * declared in the triangulation library (TriaAccessor). It enables the user - * to access the degrees of freedom on lines, quads, or hexes. The first - * template argument of this class determines the dimensionality of the object - * under consideration: 1 for lines, 2 for quads, and 3 for hexes. From the - * second template argument we can deduce the dimensionality of the - * triangulation to which this object belongs as well as the dimensionality of - * the space in which it is embedded. The second argument also denotes the - * type of DoF handler we should work on. It can either be - * ::DoFHandler or hp::DoFHandler. + * This class follows mainly the route laid out by the accessor + * library declared in the triangulation library (TriaAccessor). It + * enables the user to access the degrees of freedom on lines, quads, + * or hexes. The first template argument of this class determines the + * dimensionality of the object under consideration: 1 for lines, 2 + * for quads, and 3 for hexes. The second argument denotes the type of + * DoF handler we should work on. It can either be ::DoFHandler or + * hp::DoFHandler. From the second template argument we also deduce + * the dimension of the Triangulation this object refers to as well as + * the dimension of the space into which it is embedded. Finally, the + * template argument level_dof_access governs the + * behavior of the function dof_indices(). See the section on Generic loops below. * - * Depending on whether the structural dimension of the object - * accessed equals the dimension on which the DoF handler object - * operates, this class is derived from CellAccessor or - * TriaAccessor. This means that, for example accessors to quads - * in 2d have access to all the mesh aspects of cells, whereas - * accessors to quads in 3d can only access things that make sense for - * faces. - * - * - *

Usage

+ *

Typedefs

* * Usage is best to happen through the typedefs to the various kinds * of iterators provided by the DoFHandler and hp::DoFHandler classes, @@ -153,6 +146,22 @@ namespace internal * template interface as well as providing easier typing (much less * complicated names!). * + *

Generic loops and the third template argument

+ * + * Many loops look very similar, whether they operate on the active + * dofs of the active cells of the Triangulation or on the level dodfs + * of a single level or the whole grid hierarchy. In order to use + * polymorphism in such loops, they access degrees of freedom through + * the function dof_indices(), which changes behavior according to the + * third template argument. If the argument is false, then the active + * dofs of active cells are accessed. If it is true, the level dofs + * are used. DoFHandler has functions, for instance begin() and + * begin_mg(), which return either type or the other. Additionally, + * they can be cast into each other, in case this is needed, since + * they access the same data. + * + * It is highly recommended to use the function dof_indices() in + * generic loops in lieu of get_dof_indices() or get_mg_dof_indices(). * *

Inheritance

* @@ -169,9 +178,9 @@ namespace internal * * @ingroup dofs * @ingroup Accessors - * @author Wolfgang Bangerth, 1998, 2006, 2008 + * @author Wolfgang Bangerth, 1998, 2006, 2008, Timo Heister, Guido Kanschat, 2012 */ -template +template class DoFAccessor : public dealii::internal::DoFAccessor::Inheritance::BaseClass { public: @@ -257,9 +266,16 @@ public: * make sense, just like the * previous one. */ - template - DoFAccessor (const DoFAccessor &); + template + DoFAccessor (const DoFAccessor &); + /** + * Copy constructor allowing to + * switch level access and active + * access. + */ + template + DoFAccessor(const DoFAccessor &); /** * @} */ @@ -276,7 +292,8 @@ public: * Implement the copy operator needed * for the iterator classes. */ - void copy_from (const DoFAccessor &a); + template + void copy_from (const DoFAccessor &a); /** * Copy operator used by the @@ -287,11 +304,18 @@ public: */ void copy_from (const TriaAccessorBase &da); + /** + * Tell the caller whether + * dof_indices() accesses active + * or level dofs. + */ + static bool is_level_cell(); + /** * Return an iterator pointing to * the the parent. */ - TriaIterator > + TriaIterator > parent () const; /** @@ -305,7 +329,7 @@ public: * Return an iterator pointing to * the the @p c-th child. */ - TriaIterator > + TriaIterator > child (const unsigned int c) const; /** @@ -317,7 +341,7 @@ public: * function returns an iterator * to itself. */ - typename dealii::internal::DoFHandler::Iterators::line_iterator + typename dealii::internal::DoFHandler::Iterators::line_iterator line (const unsigned int i) const; /** @@ -329,7 +353,7 @@ public: * function returns an iterator * to itself. */ - typename dealii::internal::DoFHandler::Iterators::quad_iterator + typename dealii::internal::DoFHandler::Iterators::quad_iterator quad (const unsigned int i) const; /** @@ -417,6 +441,9 @@ public: void get_dof_indices (std::vector &dof_indices, const unsigned int fe_index = DH::default_fe_index) const; + void get_mg_dof_indices (const int level, std::vector &dof_indices, const unsigned int fe_index = DH::default_fe_index) const; + void set_mg_dof_indices (const int level, const std::vector &dof_indices, const unsigned int fe_index = DH::default_fe_index); + /** * Global DoF index of the i * degree associated with the @p vertexth @@ -457,6 +484,8 @@ public: const unsigned int i, const unsigned int fe_index = DH::default_fe_index) const; + unsigned int mg_vertex_dof_index (const int level, const unsigned int vertex, const unsigned int i, const unsigned int fe_index = DH::default_fe_index) const; + /** * Index of the ith degree * of freedom of this object. @@ -491,29 +520,31 @@ public: * freedom, and fe_index has to * match the result of * active_fe_index(). - * - * @note While the get_dof_indices() - * function returns an array that - * contains the indices of all degrees of - * freedom that somehow live on this - * object (i.e. on the vertices, edges or - * interior of this object), the current - * dof_index() function only considers - * the DoFs that really belong to this - * particular object's interior. In other - * words, as an example, if the current - * object refers to a quad (a cell in 2d, - * a face in 3d) and the finite element - * associated with it is a bilinear one, - * then the get_dof_indices() will return - * an array of size 4 while dof_index() - * will produce an exception because no - * degrees are defined in the interior of - * the face. + * + * @note While the get_dof_indices() + * function returns an array that + * contains the indices of all degrees of + * freedom that somehow live on this + * object (i.e. on the vertices, edges or + * interior of this object), the current + * dof_index() function only considers + * the DoFs that really belong to this + * particular object's interior. In other + * words, as an example, if the current + * object refers to a quad (a cell in 2d, + * a face in 3d) and the finite element + * associated with it is a bilinear one, + * then the get_dof_indices() will return + * an array of size 4 while dof_index() + * will produce an exception because no + * degrees are defined in the interior of + * the face. */ unsigned int dof_index (const unsigned int i, const unsigned int fe_index = DH::default_fe_index) const; + unsigned int mg_dof_index (const int level, const unsigned int i) const; + /** * @} */ @@ -646,17 +677,41 @@ protected: * to be accessed. */ DH *dof_handler; - - /** - * Compare for equality. - */ - bool operator == (const DoFAccessor &) const; - +public: /** - * Compare for inequality. - */ - bool operator != (const DoFAccessor &) const; - + * Compare for equality. Return + * true if the two + * accessors refer to the same object. + * + * The template parameters of this + * function allow for a comparison + * of very different + * objects. Therefore, some of them + * are disabled. + * Namely, if the dimension, or the dof + * handler of the two objects + * differ, an exception is + * generated. It can be + * expected that this is an + * unwanted comparison. + * + * The template parameter + * level_dof_access2 is ignored, such + * that an iterator with level + * access can be equal to one with + * access to the active degrees of + * freedom. + */ + template + bool operator == (const DoFAccessor &) const; + + /** + * Compare for inequality. The + * boolean not of operator==(). + */ + template + bool operator != (const DoFAccessor &) const; +protected: /** * Reset the DoF handler pointer. */ @@ -702,6 +757,8 @@ protected: const unsigned int index, const unsigned int fe_index = DH::default_fe_index) const; + void set_mg_dof_index (const int level, const unsigned int i, const unsigned int index) const; + /** * Set the global index of the i * degree on the @p vertex-th vertex of @@ -743,13 +800,15 @@ protected: const unsigned int index, const unsigned int fe_index = DH::default_fe_index) const; + void set_mg_vertex_dof_index (const int level, const unsigned int vertex, const unsigned int i, const unsigned int index, const unsigned int fe_index = DH::default_fe_index) const; + /** * Iterator classes need to be friends * because they need to access operator== * and operator!=. */ template friend class TriaRawIterator; - + template friend class DoFAccessor; private: /** @@ -769,8 +828,8 @@ private: * will give a linker error if used * anyway. */ - DoFAccessor & - operator = (const DoFAccessor &da); + DoFAccessor & + operator = (const DoFAccessor &da); /** * Make the DoFHandler class a friend so @@ -784,6 +843,7 @@ private: friend struct dealii::internal::DoFHandler::Implementation; friend struct dealii::internal::hp::DoFHandler::Implementation; friend struct dealii::internal::DoFCellAccessor::Implementation; + friend struct dealii::internal::DoFAccessor::Implementation; }; @@ -798,8 +858,8 @@ private: * * @author Wolfgang Bangerth, 2010 */ -template