]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finally add Yaqi's tutorial program, as step-28. Change references to it.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 May 2007 02:01:17 +0000 (02:01 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 May 2007 02:01:17 +0000 (02:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@14709 0785d39b-7218-0410-832d-ea1e28bc413d

23 files changed:
deal.II/base/include/base/parameter_handler.h
deal.II/deal.II/include/grid/grid_generator.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/include/grid/tria_accessor.h
deal.II/examples/step-28/Makefile [new file with mode: 0644]
deal.II/examples/step-28/doc/data-q1 [new file with mode: 0644]
deal.II/examples/step-28/doc/data-q2 [new file with mode: 0644]
deal.II/examples/step-28/doc/data-q3 [new file with mode: 0644]
deal.II/examples/step-28/doc/data-q4 [new file with mode: 0644]
deal.II/examples/step-28/doc/data-q5 [new file with mode: 0644]
deal.II/examples/step-28/doc/data-q6 [new file with mode: 0644]
deal.II/examples/step-28/doc/gnuplot.1 [new file with mode: 0644]
deal.II/examples/step-28/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-28/doc/intro.tex [new file with mode: 0644]
deal.II/examples/step-28/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-28/doc/step-28.convergence.png [new file with mode: 0755]
deal.II/examples/step-28/doc/step-28.error-vs-dofs.png [new file with mode: 0644]
deal.II/examples/step-28/doc/step-28.error-vs-time.png [new file with mode: 0644]
deal.II/examples/step-28/doc/step-28.grid-0.9.order2.png [new file with mode: 0755]
deal.II/examples/step-28/doc/step-28.grid-1.9.order2.png [new file with mode: 0755]
deal.II/examples/step-28/doc/step-28.solution-0.9.order2.png [new file with mode: 0755]
deal.II/examples/step-28/doc/step-28.solution-1.9.order2.png [new file with mode: 0755]
deal.II/examples/step-28/step-28.cc [new file with mode: 0644]

index fe205aecf5e4858073c68481442b43a02960e0d6..37eb5e4dfb4b3a40e5673fc25b3c39ac89871fbb 100644 (file)
@@ -772,7 +772,7 @@ namespace Patterns
  *
  *   The ParameterHandler class is discussed in detail in the @ref
  *   step_19 "step-19" example program, and is used in more realistic
- *   situation in @ref step_22 "step-22".
+ *   situation in @ref step_28 "step-28".
  *   
  *   <h3>Declaring entries</h3>
  *   
index bad47a009b977895527f96a6d47e9b250c29c5fb..1ac18b463324c3ef04e79da8203906f875ab8cb3 100644 (file)
@@ -268,7 +268,7 @@ class GridGenerator
                                      *
                                      * @note For an example of the
                                      * use of this function see the
-                                     * @ref step_22 "step-22"
+                                     * @ref step_28 "step-28"
                                      * tutorial program.
                                      */
     template <int dim>
index 0c92e0b52084b7e00a359da280846858a1f19226..b98c15cd1b796020bab30b69e52d681c7f9382fa 100644 (file)
@@ -849,7 +849,7 @@ namespace internal
  *   This material_id may be set upon construction of a
  *   triangulation (through the CellData data structure), or later
  *   through use of cell iterators. For a typical use of this
- *   functionality, see the @ref step_22 "step-22" tutorial
+ *   functionality, see the @ref step_28 "step-28" tutorial
  *   program. The functions of the GridGenerator namespace typically
  *   set the material ID of all cells to zero. When reading a
  *   triangulation, the material id must be specified in the input
index afe9f67b19f959488ee4d7a7b2d8ee21ea0f7a70..26bfa2b015f251c350c60689d91d76279b272c15 100644 (file)
@@ -3439,8 +3439,8 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
                                      * cell.
                                      *
                                      * For a typical use of this
-                                     * function, see the @ref step_22
-                                     * "step-22" tutorial program.
+                                     * function, see the @ref step_28
+                                     * "step-28" tutorial program.
                                      */
     unsigned char material_id () const;
 
@@ -3449,8 +3449,8 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
                                      * cell.
                                      *
                                      * For a typical use of this
-                                     * function, see the @ref step_22
-                                     * "step-22" tutorial program.
+                                     * function, see the @ref step_28
+                                     * "step-28" tutorial program.
                                      */
     void set_material_id (const unsigned char new_material_id) const;
 
diff --git a/deal.II/examples/step-28/Makefile b/deal.II/examples/step-28/Makefile
new file mode 100644 (file)
index 0000000..c620808
--- /dev/null
@@ -0,0 +1,154 @@
+# $Id: Makefile,v 1.20 2005/11/23 15:54:25 wolf Exp $
+
+
+# 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 = project
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<F)
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+       @echo ==============optimized===== $(<F)
+       @$(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: 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 to created 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/*/*.h)
+       @echo ============================ Remaking $@
+       @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
+               > $@ \
+         || (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-28/doc/data-q1 b/deal.II/examples/step-28/doc/data-q1
new file mode 100644 (file)
index 0000000..c06f6d7
--- /dev/null
@@ -0,0 +1,10 @@
+   0 2450  0.907410259645 4.077380000000
+   1 4289  0.907042395776 10.572393000000
+   2 8225  0.906937183449 25.223165000000
+   3 14757  0.906900660844 53.679839000000
+   4 26236  0.906886610239 108.345529000000
+   5 62804  0.906857043552 253.919398000000
+   6 86389  0.906844968762 466.115140000000
+   7 165629  0.906839125266 900.314131000000
+   8 301264  0.906836698748 1776.255968000000
+   9 444166  0.906835969493 3026.378920000000
diff --git a/deal.II/examples/step-28/doc/data-q2 b/deal.II/examples/step-28/doc/data-q2
new file mode 100644 (file)
index 0000000..df6c0e0
--- /dev/null
@@ -0,0 +1,10 @@
+   0 9522  0.906841960371 17.634319000000
+   1 15428  0.906837901031 33.703876000000
+   2 23566  0.906836075928 65.945975000000
+   3 33930  0.906835500110 120.341705000000
+   4 60358  0.906835001796 246.337551000000
+   5 86798  0.906834858174 403.059725000000
+   6 98538  0.906834824060 545.820023000000
+   7 117254  0.906834787555 734.265375000000
+   8 147302  0.906834761604 1014.816724000000
+   9 194442  0.906834746216 1366.371280000000
diff --git a/deal.II/examples/step-28/doc/data-q3 b/deal.II/examples/step-28/doc/data-q3
new file mode 100644 (file)
index 0000000..4b9d506
--- /dev/null
@@ -0,0 +1,10 @@
+   0 21218  0.906835491999 91.037161000000
+   1 34067  0.906834908620 131.201055000000
+   2 35069  0.906834859254 187.733460000000
+   3 50476  0.906834817384 279.669484000000
+   4 73375  0.906834783106 408.506898000000
+   5 90529  0.906834744225 590.097292000000
+   6 114417  0.906834730885 852.591387000000
+   7 151059  0.906834726355 1161.735390000000
+   8 215292  0.906834723745 1576.734300000000
+   9 275614  0.906834722268 2160.737518000000
diff --git a/deal.II/examples/step-28/doc/data-q4 b/deal.II/examples/step-28/doc/data-q4
new file mode 100644 (file)
index 0000000..15a4abe
--- /dev/null
@@ -0,0 +1,10 @@
+   0 37538  0.906834753543 373.753181000000
+   1 64160  0.906834728274 517.920264000000
+   2 88822  0.906834726164 1138.467927000000
+   3 157266  0.906834722755 2460.279981000000
+   4 175774  0.906834722205 4364.633475000000
+   5 227672  0.906834721746 4964.747243999999
+   6 264428  0.906834721654 5461.831674999999
+   7 328106  0.906834721515 6010.156317000000
+   8 395036  0.906834721416 12474.065654000000
+   9 455498  0.906834721380 19863.514286999998
diff --git a/deal.II/examples/step-28/doc/data-q5 b/deal.II/examples/step-28/doc/data-q5
new file mode 100644 (file)
index 0000000..c4aa83f
--- /dev/null
@@ -0,0 +1,10 @@
+   0 58482  0.906834731399 1121.873449000000
+   1 97921  0.906834723230 1409.121780000000
+   2 98655  0.906834723009 2989.811479000000
+   3 131934  0.906834722626 3469.306585000000
+   4 188259  0.906834722379 4011.864103000000
+   5 255493  0.906834721986 4578.781919000000
+   6 271824  0.906834721755 10632.535608000000
+   7 307409  0.906834721567 11699.187453000000
+   8 316691  0.906834721509 12528.627359000000
+   9 363748  0.906834721439 13773.087172000000
diff --git a/deal.II/examples/step-28/doc/data-q6 b/deal.II/examples/step-28/doc/data-q6
new file mode 100644 (file)
index 0000000..4c5094d
--- /dev/null
@@ -0,0 +1,15 @@
+   0 84050  0.906834722842531 3184.047951000000012
+   1 140828  0.906834721670094 3965.101212999999916
+   2 158782  0.906834721605666 8725.235561999999845
+   3 230198  0.906834721530541 16501.528386000001774
+   4 320258  0.906834721476372 18215.300853000000643
+   5 385470  0.906834721418083 21774.706741000001784
+   6 438090  0.906834721389658 25271.607131000000663
+   7 462942  0.906834721380828 27951.595711000001756
+   8 502434  0.906834721371133 33127.243892999998934
+   9 522814  0.906834721367467 36477.579565000000002
+   10 552538  0.906834721365696 39634.395656000000599
+   11 570310  0.906834721364977 43278.281700999999885
+   12 593214  0.906834721364745 47564.244135000000824
+   13 654882  0.906834721363022 53225.507492000004277
+   14 716338  0.906834721360716 60232.646245000003546
diff --git a/deal.II/examples/step-28/doc/gnuplot.1 b/deal.II/examples/step-28/doc/gnuplot.1
new file mode 100644 (file)
index 0000000..42a47f1
--- /dev/null
@@ -0,0 +1,26 @@
+set term png
+set logsc xy
+set da sty lp
+
+set ylabel "Error"
+
+set output "error-vs-dofs.png"
+set xlabel "Degrees of freedom"
+
+pl "data-q1" us 2:(abs($3-0.906834721360716)) title "Q1 elements", \
+   "data-q2" us 2:(abs($3-0.906834721360716)) title "Q2 elements", \
+   "data-q3" us 2:(abs($3-0.906834721360716)) title "Q3 elements", \
+   "data-q4" us 2:(abs($3-0.906834721360716)) title "Q4 elements", \
+   "data-q5" us 2:(abs($3-0.906834721360716)) title "Q5 elements", \
+   "data-q6" us 2:(abs($3-0.906834721360716)) title "Q6 elements"
+
+
+set output "error-vs-time.png"
+set xlabel "Time in seconds"
+
+pl "data-q1" us 4:(abs($3-0.906834721360716)) title "Q1 elements", \
+   "data-q2" us 4:(abs($3-0.906834721360716)) title "Q2 elements", \
+   "data-q3" us 4:(abs($3-0.906834721360716)) title "Q3 elements", \
+   "data-q4" us 4:(abs($3-0.906834721360716)) title "Q4 elements", \
+   "data-q5" us 4:(abs($3-0.906834721360716)) title "Q5 elements", \
+   "data-q6" us 4:(abs($3-0.906834721360716)) title "Q6 elements"
diff --git a/deal.II/examples/step-28/doc/intro.dox b/deal.II/examples/step-28/doc/intro.dox
new file mode 100644 (file)
index 0000000..f693a95
--- /dev/null
@@ -0,0 +1 @@
+<a name="Intro"></a> <h1>Introduction</h1>
diff --git a/deal.II/examples/step-28/doc/intro.tex b/deal.II/examples/step-28/doc/intro.tex
new file mode 100644 (file)
index 0000000..4516336
--- /dev/null
@@ -0,0 +1,634 @@
+\documentclass{article}
+\usepackage{amssymb,amsmath}
+\makeatletter
+\newcommand{\rmnum}[1]{\romannumeral #1}
+\newcommand{\Rmnum}[1]{\expandafter\@slowromancap\romannumeral #1@}
+\makeatother
+\begin{document}
+
+What is new in this example:
+\begin{enumerate}
+\item Solve multigroup neutron diffusion problem with multiple different meshes
+\item Solve an eigenvalue problem
+\item Setting up complicated material properties for nuclear fuel assemblies
+\end{enumerate}
+
+\subsection{Introduction}
+
+In this example, we intend to solve the multigroup diffusion approximation of
+the neutron transport equation. Essentially, the way to view this is as follows: In a
+nuclear reactor, neutrons are speeding around at different energies, get
+absorbed or scattered, or start a new fission
+event. If viewed at long enough length scales, the movement of neutrons can be
+considered a diffusion process. 
+
+A mathematical description of this would group neutrons into energy bins, and
+consider the balance equations for the neutron fluxes in each of these
+bins, or energy groups. The scattering, absorption, and fission events would
+then be operators within the diffusion equation describing the neutron
+fluxes. Assume we have energy groups $g=1,\ldots,G$, where by convention we
+assume that the neutrons with the highest energy are in group 1 and those with
+the lowest energy in group $G$. Then the neutron flux of each group satisfies the
+following equations:
+\begin{eqnarray*}
+\frac 1{v_g}\frac{\partial \phi_g(x,t)}{\partial t}
+&=&
+\nabla \cdot(D_g(x) \nabla \phi_g(x,t))
+-
+\Sigma_{r,g}(x)\phi_g(x,t)
+\\
+&& \qquad
++
+\chi_g\sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}(x,t)
++
+\sum_{g'\ne g}\Sigma_{s,g'\to g}(x)\phi_{g'}(x,t)
++
+s_{\mathrm{ext},g}(x,t)
+\end{eqnarray*}
+augmented by appropriate boundary conditions. Here, $v_g$ is the velocity of
+neutrons within group $g$. In other words, the change in
+time in flux of neutrons in group $g$ is governed by the following
+processes:
+\begin{itemize}
+\item Diffusion $\nabla \cdot(D_g(x) \nabla \phi_g(x,t))$. Here, $D_g$ is the
+  (spatially variable) diffusion coefficient. 
+\item Absorption $\Sigma_{r,g}(x)\phi_g(x,t)$ (note the
+  negative sign). The coefficient $\Sigma_{r,g}$ is called the \textit{removal
+  cross section}.
+\item Nuclear fission $\chi_g\sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}(x,t)$. 
+  The production of neutrons of energy $g$ is
+  proportional to the flux of neutrons of energy $g'$ times the
+  probability $\Sigma_{f,g'}$ that neutrons of energy $g'$ cause a fission
+  event times the number $\nu$ of neutrons produced in each fission event
+  times the probability that a neutron produced in this event has energy
+  $g$. $\nu\Sigma_{f,g'}$ is called the \textit{fission cross section} and
+  $\chi_g$ the \textit{fission spectrum}. We will denote the term
+  $\chi_g\nu\Sigma_{f,g'}$ as the \textit{fission distribution cross
+    section} in the program.
+\item Scattering $\sum_{g'\ne g}\Sigma_{s,g'\to g}(x)\phi_{g'}(x,t)$ 
+  of neutrons of energy $g'$ producing neutrons
+  of energy $g$. $\Sigma_{s,g'\to g}$ is called the \textit{scattering cross
+    section}. The case of elastic, in-group scattering $g'=g$ exists, too, but
+  we subsume this into the removal cross section. The case $g'<g$ is called
+  down-scattering, since a neutron loses energy in such an event. On the
+  other hand, $g'>g$ corresponds to up-scattering: a neutron gains energy in
+  a scattering event from the thermal motion of the atoms surrounding it;
+  up-scattering is therefore only an important process for neutrons with
+  kinetic energies that are already on the same order as the thermal kinetic
+  energy (i.e. in the sub $eV$ range).
+\item An extraneous source $s_{\mathrm{ext},g}$.
+\end{itemize}
+
+For realistic simulations in reactor analysis, one may want to split the
+continuous spectrum of neutron energies into many energy groups, often up to
+100. However, if neutron energy spectra are known well enough for some type of
+reactor (for example Pressurized Water Reactors, PWR), it is possible to obtain
+satisfactory results with only 2 energy groups. 
+
+In the program shown in this tutorial program, we provide the structure to
+compute with as many energy groups as desired. However, to keep computing
+times moderate and in order to avoid tabulating hundreds of coefficients, we
+only provide the coefficients for above equations for a two-group simulation,
+i.e. $g=1,2$. We do, however, consider a realistic situation by assuming that
+the coefficients are not constant, but rather depend on the materials that are
+assembled into reactor fuel assemblies in rather complicated ways (see
+below). 
+
+
+\subsection{The eigenvalue problem}
+
+If we consider all energy groups at once, we may write above equations in the
+following operator form:
+\begin{equation}
+\frac 1v \frac{\partial \phi}{\partial t}
+=
+-L\phi
++
+F\phi
++
+X\phi
++
+s_{\mathrm{ext}},
+\end{equation}
+where $L,F,X$ are sinking, fission, and scattering operators,
+respectively. $L$ here includes both the diffusion and removal terms. Note
+that $L$ is symmetric, whereas $F$ and $X$ are not.
+
+It is well known that this equation admits a stable solution if all
+eigenvalues of the operator $-L+F+X$ are negative. This can be readily seen by
+multiplying the equation by $\phi$ and integrating over the domain, leading to
+\begin{equation}
+  \frac 1{2v} \frac{\partial}{\partial t}  \|\phi\|^2 = ((-L+F+X)\phi,\phi).
+\end{equation}
+Stability means that the solution does not grow, i.e. we want the left hand
+side to be less than zero, which is the case if the eigenvalues of the
+operator on the right are all negative. For obvious reasons, it is
+not very desirable if a nuclear reactor produces neutron fluxes that grow
+exponentially, so eigenvalue analyses are the bread-and-butter of nuclear
+engineers. The main point of the program is therefore to consider the
+eigenvalue problem
+\begin{equation}
+  (L-F-X) \phi = \lambda \phi,
+\end{equation}
+where we want to make sure that all eigenvalues are positive. Note that $L$,
+being the diffusion operator plus the absorption (removal), is positive
+definite; the condition that all eigenvalues are positive therefore means that
+we want to make sure that fission and inter-group scattering are weak enough
+to not shift the spectrum into the negative.
+
+In nuclear engineering, one typically looks at a slightly different
+formulation of the eigenvalue problem. To this end, we do not just multiply
+with $\phi$ and integrate, but rather multiply with $\phi(L-X)^{-1}$. We then
+get the following evolution equation:
+\begin{equation}
+  \frac 1{2v} \frac{\partial}{\partial t}  \|\phi\|^2_{(L-X)^{-1}} = ((L-X)^{-1}(-L+F+X)\phi,\phi).
+\end{equation}
+Stability is the guaranteed if the eigenvalues of the following problem are
+all negative:
+\begin{equation}
+  (L-X)^{-1}(-L+F+X)\phi = \lambda_F \phi,
+\end{equation}
+which is equivalent to the eigenvalue problem
+\begin{equation}
+  (L-X)\phi = \frac 1{\lambda_F+1} F \phi.
+\end{equation}
+The typical formulation in nuclear engineering is to write this as
+\begin{equation}
+  (L-X) \phi = \frac 1{k_{\mathrm{eff}}} F \phi,
+\end{equation}
+where $k_{\mathrm{eff}}=\frac 1{\lambda^F+1}$.
+Intuitively, $k_{\mathrm{eff}}$ is something like the multiplication
+factor for neutrons per typical time scale and should be less than or equal to
+one for stable operation of a reactor: if it is less than one, the chain reaction will
+die down, whereas nuclear bombs for example have a $k$-eigenvalue larger than
+one. A stable reactor should have $k_{\mathrm{eff}}=1$.
+
+[For those who wonder how this can be achieved in practice without
+inadvertently getting slightly larger than one and triggering a nuclear bomb:
+first, fission processes happen on different time scales. While most neutrons
+are releases very quickly after a fission event, a small number of neutrons
+are only released by daughter nuclei after several further decays, up to 10-60
+seconds after the fission was initiated. If one is therefore slightly beyond
+$k_{\mathrm{eff}}=1$, one therefore has many seconds to react until all the
+neutrons created in fission re-enter the fission cycle. Nevertheless, control
+rods in nuclear reactors absorbing neutrons -- and therefore reducing
+$k_{\mathrm{eff}}$ -- are designed in such a way that they are all the way in
+the reactor in at most 2 seconds. 
+
+One therefore has on the order of 10-60 seconds to regulate the nuclear reaction
+if $k_{\mathrm{eff}}$ should be larger than one for some time, as indicated by
+a growing neutron flux. Regulation can be achieved by continuously monitoring
+the neutron flux, and if necessary increase or reduce neutron flux by moving
+neutron-absorbing control rods a few millimeters into or out of the
+reactor. On a longer scale, the water cooling the reactor contains boron, a
+good neutron absorber. Every few hours, boron concentrations are adjusted by
+adding boron or diluting the coolant.
+
+Finally, some of the absorption and scattering reactions have some
+stability built in; for example, higher neutron fluxes result in locally
+higher temperatures, which lowers the density of water and therefore reduces
+the number of scatterers that are necessary to moderate neutrons from high to
+low energies before they can start fission events themselves.]
+
+In this tutorial program, we solve above $k$-eigenvalue problem for two energy
+groups, and we are looking for the largest multiplication factor
+$k_{\mathrm{eff}}$, which is proportional to the inverse of the minimum
+eigenvalue plus one. To solve the eigenvalue problem, we generally
+use a modified version of the \emph{inverse power method}. The algorithm looks
+like this:
+
+\begin{enumerate}
+\item Initialize $\phi_g$ and $k_{\mathrm{eff}}$ with $\phi_g^{(0)}$
+  and $k_{\mathrm{eff}}^{(0)}$ and let $n=1$. 
+
+\item Define the so-called \textit{fission source} by
+  \begin{equation}
+    s_f^{(n-1)}(x)
+    =
+    \frac{1}{k_{\mathrm{eff}}^{(n-1)}}
+    \sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}^{(n-1)}(x).
+  \end{equation}
+
+\item Solve for all group fluxes $\phi_g,g=1,\ldots,G$ using
+  \begin{equation}
+    -\nabla \cdot D_g\nabla \phi_g^{(n)}
+    +
+    \Sigma_{r,g}\phi_g^{(n)}
+    = 
+    \chi_g s_f^{(n-1)}
+    + 
+    \sum_{g'< g} \Sigma_{s,g'\to g} \phi_{g'}^{(n)}
+    +
+    \sum_{g'> g}\Sigma_{s,g'\to g}\phi_{g'}^{(n-1)}.
+  \end{equation}
+
+\item Update
+  \begin{equation}
+    k_{\mathrm{eff}}^{(n)}
+    =
+    \sum_{g'=1}^G
+    \int_{\Omega}\nu\Sigma_{f,g'}(x)
+    \phi_{g'}^{(n)}(x)dx.
+  \end{equation}
+
+\item Compare $k_{\mathrm{eff}}^{(n)}$ with $k_{\mathrm{eff}}^{(n-1)}$.
+  If the change greater than a prescribed tolerance then set $n=n+1$ repeat
+  the iteration starting at step 2, otherwise end the iteration.
+\end{enumerate}
+
+Note that in this scheme, we do not solve group fluxes exactly in each power
+iteration, but rather consider previously compute $\phi_{g'}^{(n)}$ only for
+down-scattering events $g'<g$. Up-scattering is only treated by using old
+iterators $\phi_{g'}^{(n-1)}$, in essence assuming that the scattering
+operator is triangular. This is physically motivated since up-scattering does
+not play a too important role in neutron scattering. In addition, practices
+shows that the inverse power iteration is stable even using this
+simplification.
+
+Note also that one can use lots of extrapolation techniques to accelerate the
+power iteration laid out above. However, none of these are implemented in this
+example.
+
+
+\subsection{Meshes and mesh refinement}
+
+One may wonder whether it is appropriate to solve for the solutions of the
+individual energy group equations on the same meshes. The question boils down
+to this: will $\phi_g$ and $\phi_{g'}$ have similar smoothness properties? If
+this is the case, then it is appropriate to use the same mesh for the two; a
+typical application could be chemical combustion, where typically the
+concentrations of all or most chemical species change rapidly within the flame
+front. As it turns out, and as will be apparent by looking at the
+graphs shown in the results section of this tutorial program, this isn't the
+case here, however: since the diffusion coefficient is different for different
+energy groups, fast neutrons (in bins with a small group number $g$) have a very
+smooth flux function, whereas slow neutrons (in bins with a large group
+number) are much more affected by the local material properties and have a
+correspondingly rough solution if the coefficient are rough as in the case we
+compute here. Consequently, we will want to use different meshes to compute
+each energy group.
+
+This has two implications that we will have to consider: First, we need to
+find a way to refine the meshes individually. Second, assembling the source
+terms for the inverse power iteration, where we have to integrate solution
+$\phi_{g'}^{(n)}$ defined on mesh $g'$ against the shape functions defined on
+mesh $g$, becomes a much more complicated task.
+
+
+\subsubsection{Mesh refinement}
+
+We use the usual paradigm: solve on a given mesh, then evaluate an error
+indicator for each cell of each mesh we have. Because it is so convenient, we
+again use the \emph{a posteriori} error estimator by Kelly, Gago, Zienkiewicz
+and Babuska which approximates the error per cell by integrating the jump of
+the gradient of the solution along the faces of each cell. Using this, we
+obtain indicators
+\begin{equation}
+\eta_{g,K}, \qquad g=1,2,\ldots,G,\qquad K\in{\mathbb T}_g,
+\end{equation}
+where ${\mathbb T}_g$ is the triangulation used in the solution of
+$\phi_g$. The question is what to do with this. For one, it is clear that
+refining only those cells with the highest error indicators might lead to bad
+results. To understand this, it is important to realize that $\eta_{g,K}$
+scales with the second derivative of $\phi_g$. In other words, if we have two
+energy groups $g=1,2$ whose solutions are equally smooth but where one is
+larger by a factor of 10,000, for example, then only the cells of that mesh
+will be refined, whereas the mesh for the solution of small magnitude will
+remain coarse. This is probably not what one wants, since we can consider both
+components of the solution equally important.
+
+In essence, we would therefore have to scale $\eta_{g,K}$ by an importance
+factor $z_g$ that says how important it is to resolve $\phi_g$ to any given
+accuracy. Such important factors can be computed using duality techniques
+(see, for example, the step-14 tutorial program, and the
+reference to the book by Bangerth and Rannacher cited there). We won't go
+there, however, and simply assume that all energy groups are equally
+important, and will therefore normalize the error indicators $\eta_{g,K}$ for
+group $g$ by the maximum of the solution $\phi_g$. We then refine the cells
+whose errors satisfy
+\begin{equation}
+  \frac{\eta_{g,K}}{\|\phi_g\|_\infty}
+  >
+  \alpha_1
+  \displaystyle{\max_{\substack{1\le g\le G\\K\in {\mathbb T}_g}}
+    \frac{\eta_{g,K}}{\|\phi_g\|_\infty}}
+\end{equation}
+and coarsen the cells where
+\begin{equation}
+  \frac{\eta_{g,K}}{\|\phi_g\|_\infty}
+  <
+  \alpha_2
+  \displaystyle{\max_{\substack{1\le g\le G\\K\in {\mathbb T}_g}}
+    \frac{\eta_{g,K}}{\|\phi_g\|_\infty}}.
+\end{equation}
+We chose $\alpha_1=0.3$ and $\alpha_2=0.01$ in the code. Note that this will,
+of course, lead to different meshes for the different energy groups.
+
+The strategy above essentially means the following: If for energy group $g$
+there are many cells $K\in {\mathbb T}_g$ on which the error is large, for
+example because the solution is globally very rough, then many cells will be
+above the threshold. On the other hand, if there are a few cells with large
+and many with small errors, for example because the solution is overall rather
+smooth except at a few places, then only the few cells with large errors will
+be refined. Consequently, the strategy allows for meshes that track the global
+smoothness properties of the corresponding solutions rather well.
+
+
+\subsubsection{Assembling terms on different meshes}
+
+As pointed out above, the multigroup refinement strategy results in
+different meshes for the different solutions $\phi_g$. So what's the problem?
+In essence it goes like this: in step 3 of the eigenvalue iteration, we have
+form the weak form for the equation to compute $\phi_g^{(n)}$ as usual by
+multiplication with test functions $\varphi_g^i$ defined on the mesh for
+energy group $g$; in the process, we have to
+compute the right hand side vector that contains terms of the following form:
+\begin{equation}
+  F_i = \int_\Omega f(x) \varphi_g^i(x) \phi_{g'}(x) \ dx,
+\end{equation}
+where $f(x)$ is one of the coefficient functions $\Sigma_{s,g'\to g}$ or
+$\nu\chi_g\Sigma_{f,g'}$ used in the right hand side 
+of eigenvalue equation. The difficulty now is that $\phi_{g'}$ is defined on
+the mesh for energy group $g'$, i.e. it can be expanded as
+$\phi_{g'}(x)=\sum_j\phi_{g'}^j \varphi_{g'}^j(x)$, with basis functions
+$\varphi_{g'}^j(x)$ defined on mesh $g'$. The contribution to the right hand
+side can therefore be written as
+\begin{equation}
+  F_i = \sum_j \left\{\int_\Omega f(x) \varphi_g^i(x) \varphi_{g'}^j(x) 
+  \ dx \right\} \phi_{g'}^j , 
+\end{equation}
+On the other hand, the test functions $\varphi_g^i(x)$ are defined on mesh
+$g$. This means that we can't just split the integral $\Omega$ into integrals
+over the cells of either mesh $g$ or $g'$, since the respectively other basis
+functions may not be defined on these cells.
+
+The solution to this problem lies in the fact that both the meshes for $g$ and
+$g'$ are derived by adaptive refinement from a common coarse mesh. We can
+therefore always find a set of cells, which we denote by ${\mathbb T}_g \cap
+{\mathbb T}_{g'}$, that satisfy the following conditions:
+\begin{itemize}
+\item the union of the cells covers the entire domain, and
+\item a cell $K \in {\mathbb T}_g \cap {\mathbb T}_{g'}$ is active on at least
+  one of the two meshes.
+\end{itemize}
+A way to construct this set is to take each cell of coarse mesh and do the
+following steps: (i) if the cell is active on either ${\mathbb T}_g$ or
+${\mathbb T}_{g'}$, then add this cell to the set; (ii) otherwise, i.e. if
+this cell has children on both meshes, then do step (i) for each of the
+children of this cell. In fact, deal.II has a function
+\texttt{GridTools::get\_finest\_common\_cells} that computes exactly this set
+of cells that are active on at least one of two meshes.
+
+With this, we can write above integral as follows:
+\begin{equation}
+  F_i 
+  = 
+  \sum_{K \in {\mathbb T}_g \cap {\mathbb T}_{g'}}
+  \sum_j \left\{\int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x) 
+  \ dx \right\} \phi_{g'}^j.
+\end{equation}
+ In the code, we
+compute the right hand side in the function
+\texttt{NeutronDiffusionProblem::assemble\_rhs}, where (among other things) we
+loop over the set of common most refined cells, calling the function
+\texttt{NeutronDiffusionProblem::assemble\_common\_cell} on each pair of
+these cells.
+
+By construction, there are now three cases to be considered:
+\begin{itemize}
+\item[(i)] The cell $K$ is active on both meshes, i.e. both the basis
+  functions $\varphi_g^i$ as well as $\varphi_{g'}^j$ are defined on $K$.
+\item[(ii)] The cell $K$ is active on mesh $g$, but not $g'$, i.e. the
+  $\varphi_g^i$  are defined on $K$, whereas the $\varphi_{g'}^j$ are defined
+  on children of $K$.
+\item[(iii)] The cell $K$ is active on mesh $g'$, but not $g$, with opposite
+  conclusions than in (ii).
+\end{itemize}
+
+To compute the right hand side above, we then need to have different code for
+these three cases, as follows:
+\begin{itemize}
+\item[(i)] If the cell $K$ is active on both meshes, then we can directly
+  evaluate the integral. In fact, we don't even have to bother with the basis
+  functions $\varphi_{g'}$, since all we need is the values of $\phi_{g'}$ at
+  the quadrature points. We can do this using the
+  \texttt{FEValues::get\_function\_values} function. This is done directly in
+  the \texttt{NeutronDiffusionProblem::assemble\_common\_cell} function.
+
+\item[(ii)] If the cell $K$ is active on mesh $g$, but not $g'$, then the
+  basis functions $\varphi_{g'}^j$ are only defined either on the children
+  $K_c,0\le c<2^{\texttt{dim}}$, or on children of these children if cell $K$
+  is refined more than once more on mesh $g'$.
+
+  Let us assume for a second that $K$ is only once more refined on mesh $g'$
+  than on mesh $g$. Using the fact that we use embedded finite element spaces
+  where each basis function on one mesh can be written as a linear combination
+  of basis functions on the next refined mesh, we can expand the restriction
+  of $\phi_g^i$ to child cell $K_c$ into the basis functions defined on that
+  child cell (i.e. on cells on which the basis functions $\varphi_{g'}^l$ are
+  defined):
+  \begin{equation}
+    \phi_g^i|_{K_c} = B_c^{il} \varphi_{g'}^l|_{K_c}.
+  \end{equation}
+  Here, and in the following, summation over indices appearing twice is
+  implied. The matrix $B_c$ is the matrix that interpolated data from a cell
+  to its $c$-th child.
+
+  Then we can write the contribution of cell $K$ to the right hand side
+  component $F_i$ as
+  \begin{eqnarray*}
+    F_i|_K
+    &=&
+    \left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x) 
+    \ dx \right\} \phi_{g'}^j 
+    \\
+    &=&
+    \left\{
+    \sum_{0\le c<2^{\texttt{dim}}}
+    B_c^{il} \int_{K_c} f(x) \varphi_{g'}^l(x) \varphi_{g'}^j(x) 
+    \ dx \right\} \phi_{g'}^j.
+  \end{eqnarray*}
+  In matrix notation, this can be written as
+  \begin{eqnarray*}
+    F_i|_K
+    =
+    \sum_{0\le c<2^{\texttt{dim}}}
+    F_i|_{K_c},
+    \qquad
+    \qquad
+    F_i|_{K_c} = B_c^{il} M_{K_c}^{lj}  \phi_{g'}^j
+    = (B_c M_{K_c})^{il} \phi_{g'}^j,
+  \end{eqnarray*}
+  where $M_{K_c}^{lj}=\int_{K_c} f(x) \varphi_{g'}^l(x) \varphi_{g'}^j(x)$ is
+  the weighted mass matrix on child $c$ of cell $K$.
+  
+  The next question is what happens if a child $K_c$ of $K$ is not
+  active. Then, we have to apply the process recursively, i.e. we have to
+  interpolate the basis functions $\varphi_g^i$ onto child $K_c$ of $K$, then
+  onto child $K_{cc'}$ of that cell, onto child $K_{cc'c''}$ of that one, etc,
+  until we find an active cell. We then have to sum up all the contributions
+  from all the children, grandchildren, etc, of cell $K$, with contributions
+  of the form
+  \begin{equation}
+    F_i|_{K_{cc'}} = (B_cB_{c'} M_{K_{cc'}})^{ij}  \phi_{g'}^j,
+  \end{equation}
+  or
+  \begin{equation}
+    F_i|_{K_{cc'c''}} = (B_c B_{c'} B_{c''}M_{K_{cc'c''}})^{ij}
+    \phi_{g'}^j,
+  \end{equation}
+  etc. We do this process recursively, i.e. if we sit on cell $K$ and see that
+  it has children on grid $g'$, then we call a function
+  \texttt{assemble\_case\_2} with an identity matrix; the function will
+  multiply it's argument from the left with the prolongation matrix; if the
+  cell has further children, it will call itself with this new matrix,
+  otherwise it will perform the integration.
+
+\item[(iii)] The last case is where $K$ is active on mesh $g'$ but not mesh
+  $g$. In that case, we have to express basis function $\varphi_{g'}^j$ in
+  terms of the basis functions defined on the children of cell $K$, rather
+  than $\varphi_g^i$ as before. This of course works in exactly the same
+  way. If the children of $K$ are active on mesh $g$, then
+  leading to the expression
+  \begin{eqnarray*}
+    F_i|_K
+    &=&
+    \left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x) 
+    \ dx \right\} \phi_{g'}^j 
+    \\
+    &=&
+    \left\{
+    \sum_{0\le c<2^{\texttt{dim}}}
+    \int_{K_c} f(x) \varphi_{g'}^i(x) B_c^{jl} \varphi_{g'}^l(x) 
+    \ dx \right\} \phi_{g'}^j.
+  \end{eqnarray*}
+  In matrix notation, this expression now reads as
+  \begin{eqnarray*}
+    F_i|_K
+    =
+    \sum_{0\le c<2^{\texttt{dim}}}
+    F_i|_{K_c},
+    \qquad
+    \qquad
+    F_i|_{K_c} = M_{K_c}^{il} B_c^{jl}  \phi_{g'}^j
+    =
+    (M_{K_c} B_c^T)^{ij} \phi_{g'}^j,
+  \end{eqnarray*}
+  and correspondingly for cases where cell $K$ is refined more than once on
+  mesh $g$:
+  \begin{equation}
+    F_i|_{K_{cc'}} = (M_{K_{cc'}} B_{c'}^T B_c^T)^{ij}  \phi_{g'}^j,
+  \end{equation}
+  or
+  \begin{equation}
+    F_i|_{K_{cc'c''}} = (M_{K_{cc'c''}} B_{c''}^T B_{c'}^T B_c^T)^{ij}
+    \phi_{g'}^j,
+  \end{equation}
+  etc. In other words, the process works in exactly the same way as before,
+  except that we have to take the transpose of the prolongation matrices and
+  need to multiply it to the mass matrix from the other side. 
+\end{itemize}
+
+
+The expressions for cases (ii) and (iii) can be understood as repeatedly
+interpolating either the left or right basis functions in the scalar product
+$(f \varphi_g^i, \varphi_{g'}^j)_K$ onto child cells, and then finally
+forming the inner product (the mass matrix) on the final cell. To make the
+symmetry in these cases more obvious, we can write them like this: for case
+(ii), we have
+\begin{equation}
+  F_i|_{K_{cc'\cdots c^{(k)}}}
+  = [B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}}]^{ij}
+    \phi_{g'}^j,
+\end{equation}
+whereas for case (iii) we get
+\begin{equation}
+  F_i|_{K_{cc'\cdots c^{(k)}}}
+  = [(B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}})^T]^{ij}
+    \phi_{g'}^j,
+\end{equation}
+
+
+
+\subsection{Description of the test case}
+
+A nuclear reactor core is composed of different types of assemblies. An
+assembly is essentially the smallest unit that can be moved in and out of a
+reactor, and is usually rectangular or square. However, assemblies are not
+fixed units, as they are assembled from a complex lattice of different fuel
+rods, control rods, and instrumentation elements that are held in place
+relative to each other by spacers that are permanently attached to the rods.
+To make things more complicated, there are different kinds of assemblies that
+are used at the same time in a reactor, where assemblies differ in the type
+and arrangement of rods they are made up of.
+
+Obviously, the arrangement of assemblies as well as the arrangement of rods
+inside them affect the distribution of neutron fluxes in the reactor (a fact
+that will be obvious by looking at the solution shown below in the results
+sections of this program). Fuel rods, for example, differ from each other in
+the enrichment of U-235 or Pu-239. Control rods, on the other hand, have zero
+fission, but nonzero scattering and absorption cross sections.
+
+This whole arrangement would make the description or spatially dependent
+material parameters very complicated. It will not become much simpler, but we
+will make one approximation: we merge the volume inhabited by each cylindrical
+rod and the surrounding water into volumes of quadratic cross section into
+so-called ``pin cells'' for which homogenized material data are obtained with
+nuclear database and knowledge of neutron spectrum. The homogenization makes
+all material data piecewise constant on the solution domain for a reactor with
+fresh fuel. Spatially dependent material parameters are then looked up for the
+quadratic assembly in which a point is located, and then for the quadratic pin
+cell within this assembly.
+
+In this tutorial program, we simulate a quarter of a reactor consisting of $4
+\times 4$ assemblies. We use symmetry (Neumann) boundary conditions to reduce
+the problem to one quarter of the domain, and consequently only simulate a
+$2\times 2$ set of assemblies. Two of them will be UO${}_2$ fuel, the other
+two of them MOX fuel. Each of these assemblies consists of $17\times 17$ rods
+of different compositions. In total, we therefore create a $34\times 34$
+lattice of rods. To make things simpler later on, we reflect this fact by
+creating a coarse mesh of $34\times 34$ cells (even though the domain is a
+square, for which we would usually use a single cell). In deal.II, each cell
+has a \texttt{material\_id} which one may use to associated each cell with a
+particular number identifying the material from which this cell's volume is
+made of; we will use this material ID to identify which of the 8 different
+kinds of rods that are used in this testcase make up a particular cell. Note
+that upon mesh refinement, the children of a cell inherit the material ID,
+making it simple to track the material even after mesh refinement.
+
+The arrangement of the rods will be clearly visible in the images shown in
+the results section. The cross sections for materials and for both energy
+groups are taken from a OECD/NEA benchmark problem. The detailed configuration
+and material data is given in the code.
+
+
+\subsection{What the program does (and how it does that)}
+
+As a coarse overview of what exactly the program does, here is the basic
+layout: starting on a coarse mesh that is the same for each energy group, we
+compute inverse eigenvalue iterations to compute the $k$-eigenvalue on a given
+set of meshes. We stop these iterations when the change in the eigenvalue
+drops below a certain tolerance, and then write out the meshes and solutions
+for each energy group for inspection by a graphics program. Because the meshes
+for the solutions are different, we have to generate a separate output file
+for each energy group, rather than being able to add all energy group
+solutions into the same file.
+
+After this, we evaluate the error indicators as explained in one of the sections
+above for each of the meshes, and refine and coarsen the cells of each mesh
+independently. Since the eigenvalue iterations are fairly expensive, we don't
+want to start all over on the new mesh; rather, we use the SolutionTransfer
+class to interpolate the solution on the previous mesh to the next one upon
+mesh refinement. A simple experiment will convince you that this is a lot
+cheaper than if we omitted this step. After doing so, we resume our eigenvalue
+iterations on the next set of meshes.
+
+The program is controlled by a parameter file, using the ParameterHandler
+class already mentioned in the step-19 example program. We will show a
+parameter file in the results section of this section. For the moment suffice
+it to say that it controls the polynomial degree of the finite elements used,
+the number of energy groups (even though all that is presently implemented are
+the coefficients for a 2-group problem), the tolerance where to stop the
+inverse eigenvalue iteration, and the number of refinement cycles we will do.
+
+\end{document}
diff --git a/deal.II/examples/step-28/doc/results.dox b/deal.II/examples/step-28/doc/results.dox
new file mode 100644 (file)
index 0000000..f7a468d
--- /dev/null
@@ -0,0 +1,500 @@
+<a name="Results"></a> <h1>Results</h1>
+
+
+The output of this program consist of the console output, a file 
+named ``convergence_table'' to record main results of mesh iteration, the eps
+files including the grids, and the solutions given in gnuplot format.
+
+When we set Polynomial_Order to 2, we got following console output:
+@code
+Cycle 0:
+   Numbers of active cells:       1156 1156 
+   Numbers of degrees of freedom: 4761 4761 
+Iter number:1 k_eff=319.375676634307 flux ratio=6.836246075631 max_thermal=1.433899030144
+Iter number:2 k_eff=0.834072546055 flux ratio=5.204601882141 max_thermal=0.004630925876
+Iter number:3 k_eff=0.862826188043 flux ratio=4.645051765984 max_thermal=0.005380396338
+Iter number:4 k_eff=0.877887920967 flux ratio=4.318030683875 max_thermal=0.006005512201
+Iter number:5 k_eff=0.887161559547 flux ratio=4.256596788174 max_thermal=0.006639443035
+Iter number:6 k_eff=0.893254525197 flux ratio=4.296498905676 max_thermal=0.007161016401
+Iter number:7 k_eff=0.897386466621 flux ratio=4.323736110066 max_thermal=0.007541125053
+Iter number:8 k_eff=0.900235644733 flux ratio=4.342491852394 max_thermal=0.007813654241
+Iter number:9 k_eff=0.902217719823 flux ratio=4.355367629620 max_thermal=0.008007335384
+Iter number:10 k_eff=0.903602785157 flux ratio=4.364212965582 max_thermal=0.008144201718
+Iter number:11 k_eff=0.904572678811 flux ratio=4.370302672219 max_thermal=0.008240563668
+Iter number:12 k_eff=0.905252379018 flux ratio=4.374506018233 max_thermal=0.008308245052
+Iter number:13 k_eff=0.905728767660 flux ratio=4.377414535866 max_thermal=0.008355707257
+Iter number:14 k_eff=0.906062594755 flux ratio=4.379431495993 max_thermal=0.008388956163
+Iter number:15 k_eff=0.906296449777 flux ratio=4.380832749068 max_thermal=0.008412232394
+Iter number:16 k_eff=0.906460217413 flux ratio=4.381807689696 max_thermal=0.008428519956
+Iter number:17 k_eff=0.906574868174 flux ratio=4.382486812297 max_thermal=0.008439913928
+Iter number:18 k_eff=0.906655112098 flux ratio=4.382960306180 max_thermal=0.008447883119
+Iter number:19 k_eff=0.906711262646 flux ratio=4.383290664907 max_thermal=0.008453456292
+Iter number:20 k_eff=0.906750547146 flux ratio=4.383521281424 max_thermal=0.008457353551
+Iter number:21 k_eff=0.906778027984 flux ratio=4.383682335660 max_thermal=0.008460078740
+Iter number:22 k_eff=0.906797249757 flux ratio=4.383794844804 max_thermal=0.008461984300
+Iter number:23 k_eff=0.906810693563 flux ratio=4.383873459628 max_thermal=0.008463316724
+Iter number:24 k_eff=0.906820095658 flux ratio=4.383928400683 max_thermal=0.008464248386
+Iter number:25 k_eff=0.906826670830 flux ratio=4.383966802041 max_thermal=0.008464899825
+Iter number:26 k_eff=0.906831268888 flux ratio=4.383993645532 max_thermal=0.008465355326
+Iter number:27 k_eff=0.906834484255 flux ratio=4.384012411177 max_thermal=0.008465673822
+Iter number:28 k_eff=0.906836732678 flux ratio=4.384025530521 max_thermal=0.008465896521
+Iter number:29 k_eff=0.906838304919 flux ratio=4.384034702833 max_thermal=0.008466052239
+Iter number:30 k_eff=0.906839404318 flux ratio=4.384041115801 max_thermal=0.008466161120
+Iter number:31 k_eff=0.906840173074 flux ratio=4.384045599636 max_thermal=0.008466237253
+Iter number:32 k_eff=0.906840710623 flux ratio=4.384048734710 max_thermal=0.008466290487
+Iter number:33 k_eff=0.906841086501 flux ratio=4.384050926767 max_thermal=0.008466327710
+Iter number:34 k_eff=0.906841349329 flux ratio=4.384052459477 max_thermal=0.008466353737
+Iter number:35 k_eff=0.906841533109 flux ratio=4.384053531173 max_thermal=0.008466371936
+Iter number:36 k_eff=0.906841661615 flux ratio=4.384054280525 max_thermal=0.008466384662
+Iter number:37 k_eff=0.906841751471 flux ratio=4.384054804489 max_thermal=0.008466393560
+Iter number:38 k_eff=0.906841814301 flux ratio=4.384055170858 max_thermal=0.008466399781
+Iter number:39 k_eff=0.906841858234 flux ratio=4.384055427034 max_thermal=0.008466404132
+Iter number:40 k_eff=0.906841888954 flux ratio=4.384055606159 max_thermal=0.008466407174
+Iter number:41 k_eff=0.906841910434 flux ratio=4.384055731409 max_thermal=0.008466409301
+Iter number:42 k_eff=0.906841925454 flux ratio=4.384055818987 max_thermal=0.008466410788
+Iter number:43 k_eff=0.906841935956 flux ratio=4.384055880225 max_thermal=0.008466411828
+Iter number:44 k_eff=0.906841943300 flux ratio=4.384055923044 max_thermal=0.008466412555
+Iter number:45 k_eff=0.906841948435 flux ratio=4.384055952984 max_thermal=0.008466413064
+Iter number:46 k_eff=0.906841952025 flux ratio=4.384055973920 max_thermal=0.008466413419
+Iter number:47 k_eff=0.906841954536 flux ratio=4.384055988559 max_thermal=0.008466413668
+Iter number:48 k_eff=0.906841956291 flux ratio=4.384055998794 max_thermal=0.008466413842
+Iter number:49 k_eff=0.906841957518 flux ratio=4.384056005952 max_thermal=0.008466413963
+Iter number:50 k_eff=0.906841958377 flux ratio=4.384056010956 max_thermal=0.008466414048
+Iter number:51 k_eff=0.906841958977 flux ratio=4.384056014456 max_thermal=0.008466414108
+Iter number:52 k_eff=0.906841959397 flux ratio=4.384056016902 max_thermal=0.008466414149
+Iter number:53 k_eff=0.906841959690 flux ratio=4.384056018613 max_thermal=0.008466414178
+Iter number:54 k_eff=0.906841959895 flux ratio=4.384056019810 max_thermal=0.008466414199
+Iter number:55 k_eff=0.906841960039 flux ratio=4.384056020646 max_thermal=0.008466414213
+Iter number:56 k_eff=0.906841960139 flux ratio=4.384056021231 max_thermal=0.008466414223
+Iter number:57 k_eff=0.906841960209 flux ratio=4.384056021640 max_thermal=0.008466414230
+Iter number:58 k_eff=0.906841960258 flux ratio=4.384056021926 max_thermal=0.008466414235
+Iter number:59 k_eff=0.906841960292 flux ratio=4.384056022126 max_thermal=0.008466414238
+Iter number:60 k_eff=0.906841960316 flux ratio=4.384056022266 max_thermal=0.008466414240
+Iter number:61 k_eff=0.906841960333 flux ratio=4.384056022364 max_thermal=0.008466414242
+Iter number:62 k_eff=0.906841960345 flux ratio=4.384056022432 max_thermal=0.008466414243
+Iter number:63 k_eff=0.906841960353 flux ratio=4.384056022480 max_thermal=0.008466414244
+Iter number:64 k_eff=0.906841960359 flux ratio=4.384056022513 max_thermal=0.008466414245
+Iter number:65 k_eff=0.906841960363 flux ratio=4.384056022537 max_thermal=0.008466414245
+Iter number:66 k_eff=0.906841960366 flux ratio=4.384056022553 max_thermal=0.008466414245
+Iter number:67 k_eff=0.906841960368 flux ratio=4.384056022564 max_thermal=0.008466414246
+Iter number:68 k_eff=0.906841960369 flux ratio=4.384056022572 max_thermal=0.008466414246
+Iter number:69 k_eff=0.906841960370 flux ratio=4.384056022578 max_thermal=0.008466414246
+Iter number:70 k_eff=0.906841960371 flux ratio=4.384056022582 max_thermal=0.008466414246
+Cycle 1:
+   Numbers of active cells:       1156 2380 
+   Numbers of degrees of freedom: 4761 10667 
+Iter number:1 k_eff=0.906838267472 flux ratio=4.385474405124 max_thermal=0.008463675976
+Iter number:2 k_eff=0.906837892433 flux ratio=4.385486158840 max_thermal=0.008463675386
+Iter number:3 k_eff=0.906837848258 flux ratio=4.385487761080 max_thermal=0.008463681343
+Iter number:4 k_eff=0.906837849549 flux ratio=4.385488316617 max_thermal=0.008463685560
+Iter number:5 k_eff=0.906837859133 flux ratio=4.385488608157 max_thermal=0.008463688398
+Iter number:6 k_eff=0.906837869078 flux ratio=4.385488780872 max_thermal=0.008463690308
+Iter number:7 k_eff=0.906837877437 flux ratio=4.385488887750 max_thermal=0.008463691599
+Iter number:8 k_eff=0.906837883930 flux ratio=4.385488955514 max_thermal=0.008463692477
+Iter number:9 k_eff=0.906837888778 flux ratio=4.385488999282 max_thermal=0.008463693077
+Iter number:10 k_eff=0.906837892317 flux ratio=4.385489027995 max_thermal=0.008463693490
+Iter number:11 k_eff=0.906837894865 flux ratio=4.385489047078 max_thermal=0.008463693775
+Iter number:12 k_eff=0.906837896682 flux ratio=4.385489059900 max_thermal=0.008463693972
+Iter number:13 k_eff=0.906837897972 flux ratio=4.385489068591 max_thermal=0.008463694108
+Iter number:14 k_eff=0.906837898882 flux ratio=4.385489074523 max_thermal=0.008463694203
+Iter number:15 k_eff=0.906837899524 flux ratio=4.385489078594 max_thermal=0.008463694269
+Iter number:16 k_eff=0.906837899975 flux ratio=4.385489081400 max_thermal=0.008463694315
+Iter number:17 k_eff=0.906837900292 flux ratio=4.385489083340 max_thermal=0.008463694347
+Iter number:18 k_eff=0.906837900514 flux ratio=4.385489084686 max_thermal=0.008463694369
+Iter number:19 k_eff=0.906837900670 flux ratio=4.385489085620 max_thermal=0.008463694385
+Iter number:20 k_eff=0.906837900779 flux ratio=4.385489086270 max_thermal=0.008463694396
+Iter number:21 k_eff=0.906837900855 flux ratio=4.385489086723 max_thermal=0.008463694404
+Iter number:22 k_eff=0.906837900909 flux ratio=4.385489087039 max_thermal=0.008463694409
+Iter number:23 k_eff=0.906837900946 flux ratio=4.385489087259 max_thermal=0.008463694413
+Iter number:24 k_eff=0.906837900972 flux ratio=4.385489087413 max_thermal=0.008463694415
+Iter number:25 k_eff=0.906837900990 flux ratio=4.385489087521 max_thermal=0.008463694417
+Iter number:26 k_eff=0.906837901003 flux ratio=4.385489087596 max_thermal=0.008463694418
+Iter number:27 k_eff=0.906837901012 flux ratio=4.385489087648 max_thermal=0.008463694419
+Iter number:28 k_eff=0.906837901018 flux ratio=4.385489087685 max_thermal=0.008463694420
+Iter number:29 k_eff=0.906837901023 flux ratio=4.385489087710 max_thermal=0.008463694420
+Iter number:30 k_eff=0.906837901026 flux ratio=4.385489087728 max_thermal=0.008463694421
+Iter number:31 k_eff=0.906837901028 flux ratio=4.385489087741 max_thermal=0.008463694421
+Iter number:32 k_eff=0.906837901030 flux ratio=4.385489087749 max_thermal=0.008463694421
+Iter number:33 k_eff=0.906837901031 flux ratio=4.385489087755 max_thermal=0.008463694421
+Iter number:34 k_eff=0.906837901031 flux ratio=4.385489087760 max_thermal=0.008463694421
+Cycle 2:
+   Numbers of active cells:       1156 4219 
+   Numbers of degrees of freedom: 4761 18805 
+Iter number:1 k_eff=0.906836032131 flux ratio=4.385463219198 max_thermal=0.008463744346
+Iter number:2 k_eff=0.906835885276 flux ratio=4.385464355771 max_thermal=0.008463756043
+Iter number:3 k_eff=0.906835925468 flux ratio=4.385465086805 max_thermal=0.008463765603
+Iter number:4 k_eff=0.906835957285 flux ratio=4.385465633530 max_thermal=0.008463771996
+Iter number:5 k_eff=0.906835986095 flux ratio=4.385465985078 max_thermal=0.008463776366
+Iter number:6 k_eff=0.906836010004 flux ratio=4.385466206119 max_thermal=0.008463779390
+Iter number:7 k_eff=0.906836028489 flux ratio=4.385466348192 max_thermal=0.008463781494
+Iter number:8 k_eff=0.906836042190 flux ratio=4.385466442031 max_thermal=0.008463782961
+Iter number:9 k_eff=0.906836052101 flux ratio=4.385466505347 max_thermal=0.008463783987
+Iter number:10 k_eff=0.906836059171 flux ratio=4.385466548687 max_thermal=0.008463784704
+Iter number:11 k_eff=0.906836064172 flux ratio=4.385466578622 max_thermal=0.008463785205
+Iter number:12 k_eff=0.906836067692 flux ratio=4.385466599413 max_thermal=0.008463785556
+Iter number:13 k_eff=0.906836070164 flux ratio=4.385466613899 max_thermal=0.008463785801
+Iter number:14 k_eff=0.906836071896 flux ratio=4.385466624011 max_thermal=0.008463785972
+Iter number:15 k_eff=0.906836073108 flux ratio=4.385466631076 max_thermal=0.008463786092
+Iter number:16 k_eff=0.906836073957 flux ratio=4.385466636016 max_thermal=0.008463786176
+Iter number:17 k_eff=0.906836074550 flux ratio=4.385466639471 max_thermal=0.008463786235
+Iter number:18 k_eff=0.906836074965 flux ratio=4.385466641887 max_thermal=0.008463786276
+Iter number:19 k_eff=0.906836075255 flux ratio=4.385466643577 max_thermal=0.008463786305
+Iter number:20 k_eff=0.906836075458 flux ratio=4.385466644759 max_thermal=0.008463786325
+Iter number:21 k_eff=0.906836075600 flux ratio=4.385466645585 max_thermal=0.008463786339
+Iter number:22 k_eff=0.906836075699 flux ratio=4.385466646163 max_thermal=0.008463786349
+Iter number:23 k_eff=0.906836075769 flux ratio=4.385466646568 max_thermal=0.008463786356
+Iter number:24 k_eff=0.906836075817 flux ratio=4.385466646850 max_thermal=0.008463786360
+Iter number:25 k_eff=0.906836075851 flux ratio=4.385466647048 max_thermal=0.008463786364
+Iter number:26 k_eff=0.906836075875 flux ratio=4.385466647186 max_thermal=0.008463786366
+Iter number:27 k_eff=0.906836075891 flux ratio=4.385466647283 max_thermal=0.008463786368
+Iter number:28 k_eff=0.906836075903 flux ratio=4.385466647351 max_thermal=0.008463786369
+Iter number:29 k_eff=0.906836075911 flux ratio=4.385466647398 max_thermal=0.008463786370
+Iter number:30 k_eff=0.906836075917 flux ratio=4.385466647431 max_thermal=0.008463786370
+Iter number:31 k_eff=0.906836075921 flux ratio=4.385466647454 max_thermal=0.008463786371
+Iter number:32 k_eff=0.906836075924 flux ratio=4.385466647470 max_thermal=0.008463786371
+Iter number:33 k_eff=0.906836075926 flux ratio=4.385466647482 max_thermal=0.008463786371
+Iter number:34 k_eff=0.906836075927 flux ratio=4.385466647489 max_thermal=0.008463786371
+Iter number:35 k_eff=0.906836075928 flux ratio=4.385466647495 max_thermal=0.008463786371
+Iter number:36 k_eff=0.906836075928 flux ratio=4.385466647499 max_thermal=0.008463786371
+Cycle 3:
+   Numbers of active cells:       1507 6133 
+   Numbers of degrees of freedom: 6629 27301 
+Iter number:1 k_eff=0.906835269231 flux ratio=4.385403797661 max_thermal=0.008463904814
+Iter number:2 k_eff=0.906835336040 flux ratio=4.385404137890 max_thermal=0.008463908915
+Iter number:3 k_eff=0.906835394509 flux ratio=4.385404254946 max_thermal=0.008463912078
+Iter number:4 k_eff=0.906835430587 flux ratio=4.385404329325 max_thermal=0.008463914341
+Iter number:5 k_eff=0.906835453664 flux ratio=4.385404388953 max_thermal=0.008463915964
+Iter number:6 k_eff=0.906835468672 flux ratio=4.385404437567 max_thermal=0.008463917123
+Iter number:7 k_eff=0.906835478615 flux ratio=4.385404475817 max_thermal=0.008463917947
+Iter number:8 k_eff=0.906835485309 flux ratio=4.385404504893 max_thermal=0.008463918530
+Iter number:9 k_eff=0.906835489870 flux ratio=4.385404526447 max_thermal=0.008463918942
+Iter number:10 k_eff=0.906835493002 flux ratio=4.385404542148 max_thermal=0.008463919232
+Iter number:11 k_eff=0.906835495166 flux ratio=4.385404553451 max_thermal=0.008463919436
+Iter number:12 k_eff=0.906835496666 flux ratio=4.385404561520 max_thermal=0.008463919579
+Iter number:13 k_eff=0.906835497709 flux ratio=4.385404567248 max_thermal=0.008463919680
+Iter number:14 k_eff=0.906835498435 flux ratio=4.385404571298 max_thermal=0.008463919750
+Iter number:15 k_eff=0.906835498941 flux ratio=4.385404574153 max_thermal=0.008463919800
+Iter number:16 k_eff=0.906835499294 flux ratio=4.385404576162 max_thermal=0.008463919834
+Iter number:17 k_eff=0.906835499541 flux ratio=4.385404577572 max_thermal=0.008463919858
+Iter number:18 k_eff=0.906835499713 flux ratio=4.385404578562 max_thermal=0.008463919875
+Iter number:19 k_eff=0.906835499833 flux ratio=4.385404579256 max_thermal=0.008463919887
+Iter number:20 k_eff=0.906835499917 flux ratio=4.385404579742 max_thermal=0.008463919895
+Iter number:21 k_eff=0.906835499976 flux ratio=4.385404580083 max_thermal=0.008463919901
+Iter number:22 k_eff=0.906835500017 flux ratio=4.385404580321 max_thermal=0.008463919905
+Iter number:23 k_eff=0.906835500046 flux ratio=4.385404580488 max_thermal=0.008463919908
+Iter number:24 k_eff=0.906835500066 flux ratio=4.385404580604 max_thermal=0.008463919910
+Iter number:25 k_eff=0.906835500080 flux ratio=4.385404580686 max_thermal=0.008463919911
+Iter number:26 k_eff=0.906835500090 flux ratio=4.385404580743 max_thermal=0.008463919912
+Iter number:27 k_eff=0.906835500097 flux ratio=4.385404580783 max_thermal=0.008463919913
+Iter number:28 k_eff=0.906835500101 flux ratio=4.385404580811 max_thermal=0.008463919914
+Iter number:29 k_eff=0.906835500105 flux ratio=4.385404580831 max_thermal=0.008463919914
+Iter number:30 k_eff=0.906835500107 flux ratio=4.385404580844 max_thermal=0.008463919914
+Iter number:31 k_eff=0.906835500109 flux ratio=4.385404580854 max_thermal=0.008463919914
+Iter number:32 k_eff=0.906835500110 flux ratio=4.385404580860 max_thermal=0.008463919914
+Iter number:33 k_eff=0.906835500111 flux ratio=4.385404580865 max_thermal=0.008463919914
+Cycle 4:
+   Numbers of active cells:       2734 10732 
+   Numbers of degrees of freedom: 12263 48095 
+Iter number:1 k_eff=0.906834846364 flux ratio=4.385381150927 max_thermal=0.008463963132
+Iter number:2 k_eff=0.906834885938 flux ratio=4.385381403919 max_thermal=0.008463966541
+Iter number:3 k_eff=0.906834926531 flux ratio=4.385381460228 max_thermal=0.008463969391
+Iter number:4 k_eff=0.906834951140 flux ratio=4.385381536948 max_thermal=0.008463971436
+Iter number:5 k_eff=0.906834967134 flux ratio=4.385381607000 max_thermal=0.008463972881
+Iter number:6 k_eff=0.906834977868 flux ratio=4.385381661804 max_thermal=0.008463973895
+Iter number:7 k_eff=0.906834985190 flux ratio=4.385381702101 max_thermal=0.008463974604
+Iter number:8 k_eff=0.906834990235 flux ratio=4.385381730930 max_thermal=0.008463975099
+Iter number:9 k_eff=0.906834993731 flux ratio=4.385381751290 max_thermal=0.008463975445
+Iter number:10 k_eff=0.906834996164 flux ratio=4.385381765579 max_thermal=0.008463975687
+Iter number:11 k_eff=0.906834997860 flux ratio=4.385381775578 max_thermal=0.008463975856
+Iter number:12 k_eff=0.906834999045 flux ratio=4.385381782564 max_thermal=0.008463975973
+Iter number:13 k_eff=0.906834999873 flux ratio=4.385381787442 max_thermal=0.008463976056
+Iter number:14 k_eff=0.906835000452 flux ratio=4.385381790849 max_thermal=0.008463976113
+Iter number:15 k_eff=0.906835000856 flux ratio=4.385381793228 max_thermal=0.008463976153
+Iter number:16 k_eff=0.906835001140 flux ratio=4.385381794889 max_thermal=0.008463976182
+Iter number:17 k_eff=0.906835001338 flux ratio=4.385381796049 max_thermal=0.008463976201
+Iter number:18 k_eff=0.906835001476 flux ratio=4.385381796860 max_thermal=0.008463976215
+Iter number:19 k_eff=0.906835001573 flux ratio=4.385381797427 max_thermal=0.008463976224
+Iter number:20 k_eff=0.906835001641 flux ratio=4.385381797823 max_thermal=0.008463976231
+Iter number:21 k_eff=0.906835001688 flux ratio=4.385381798099 max_thermal=0.008463976236
+Iter number:22 k_eff=0.906835001721 flux ratio=4.385381798293 max_thermal=0.008463976239
+Iter number:23 k_eff=0.906835001744 flux ratio=4.385381798428 max_thermal=0.008463976241
+Iter number:24 k_eff=0.906835001760 flux ratio=4.385381798523 max_thermal=0.008463976243
+Iter number:25 k_eff=0.906835001772 flux ratio=4.385381798589 max_thermal=0.008463976244
+Iter number:26 k_eff=0.906835001780 flux ratio=4.385381798635 max_thermal=0.008463976245
+Iter number:27 k_eff=0.906835001785 flux ratio=4.385381798667 max_thermal=0.008463976246
+Iter number:28 k_eff=0.906835001789 flux ratio=4.385381798690 max_thermal=0.008463976246
+Iter number:29 k_eff=0.906835001792 flux ratio=4.385381798706 max_thermal=0.008463976246
+Iter number:30 k_eff=0.906835001794 flux ratio=4.385381798717 max_thermal=0.008463976246
+Iter number:31 k_eff=0.906835001795 flux ratio=4.385381798724 max_thermal=0.008463976246
+Iter number:32 k_eff=0.906835001796 flux ratio=4.385381798730 max_thermal=0.008463976247
+Iter number:33 k_eff=0.906835001796 flux ratio=4.385381798734 max_thermal=0.008463976247
+Cycle 5:
+   Numbers of active cells:       3928 15598 
+   Numbers of degrees of freedom: 17501 69297 
+Iter number:1 k_eff=0.906834756419 flux ratio=4.384853631027 max_thermal=0.008464995625
+Iter number:2 k_eff=0.906834805316 flux ratio=4.384853670539 max_thermal=0.008464997340
+Iter number:3 k_eff=0.906834826349 flux ratio=4.384853696977 max_thermal=0.008464998471
+Iter number:4 k_eff=0.906834837235 flux ratio=4.384853725424 max_thermal=0.008464999257
+Iter number:5 k_eff=0.906834843989 flux ratio=4.384853750658 max_thermal=0.008464999812
+Iter number:6 k_eff=0.906834848438 flux ratio=4.384853770680 max_thermal=0.008465000204
+Iter number:7 k_eff=0.906834851444 flux ratio=4.384853785717 max_thermal=0.008465000479
+Iter number:8 k_eff=0.906834853502 flux ratio=4.384853796687 max_thermal=0.008465000673
+Iter number:9 k_eff=0.906834854923 flux ratio=4.384853804560 max_thermal=0.008465000809
+Iter number:10 k_eff=0.906834855908 flux ratio=4.384853810156 max_thermal=0.008465000904
+Iter number:11 k_eff=0.906834856593 flux ratio=4.384853814111 max_thermal=0.008465000971
+Iter number:12 k_eff=0.906834857071 flux ratio=4.384853816895 max_thermal=0.008465001018
+Iter number:13 k_eff=0.906834857404 flux ratio=4.384853818850 max_thermal=0.008465001050
+Iter number:14 k_eff=0.906834857636 flux ratio=4.384853820222 max_thermal=0.008465001073
+Iter number:15 k_eff=0.906834857799 flux ratio=4.384853821183 max_thermal=0.008465001089
+Iter number:16 k_eff=0.906834857912 flux ratio=4.384853821856 max_thermal=0.008465001100
+Iter number:17 k_eff=0.906834857991 flux ratio=4.384853822327 max_thermal=0.008465001108
+Iter number:18 k_eff=0.906834858047 flux ratio=4.384853822657 max_thermal=0.008465001114
+Iter number:19 k_eff=0.906834858085 flux ratio=4.384853822888 max_thermal=0.008465001117
+Iter number:20 k_eff=0.906834858113 flux ratio=4.384853823049 max_thermal=0.008465001120
+Iter number:21 k_eff=0.906834858132 flux ratio=4.384853823162 max_thermal=0.008465001122
+Iter number:22 k_eff=0.906834858145 flux ratio=4.384853823241 max_thermal=0.008465001123
+Iter number:23 k_eff=0.906834858154 flux ratio=4.384853823296 max_thermal=0.008465001124
+Iter number:24 k_eff=0.906834858161 flux ratio=4.384853823335 max_thermal=0.008465001125
+Iter number:25 k_eff=0.906834858165 flux ratio=4.384853823362 max_thermal=0.008465001125
+Iter number:26 k_eff=0.906834858168 flux ratio=4.384853823380 max_thermal=0.008465001126
+Iter number:27 k_eff=0.906834858170 flux ratio=4.384853823394 max_thermal=0.008465001126
+Iter number:28 k_eff=0.906834858172 flux ratio=4.384853823403 max_thermal=0.008465001126
+Iter number:29 k_eff=0.906834858173 flux ratio=4.384853823409 max_thermal=0.008465001126
+Iter number:30 k_eff=0.906834858174 flux ratio=4.384853823414 max_thermal=0.008465001126
+Cycle 6:
+   Numbers of active cells:       4486 17755 
+   Numbers of degrees of freedom: 19933 78605 
+Iter number:1 k_eff=0.906834797396 flux ratio=4.384850615384 max_thermal=0.008465007174
+Iter number:2 k_eff=0.906834814661 flux ratio=4.384850631692 max_thermal=0.008465007547
+Iter number:3 k_eff=0.906834818206 flux ratio=4.384850637654 max_thermal=0.008465007740
+Iter number:4 k_eff=0.906834820171 flux ratio=4.384850642364 max_thermal=0.008465007876
+Iter number:5 k_eff=0.906834821415 flux ratio=4.384850646407 max_thermal=0.008465007974
+Iter number:6 k_eff=0.906834822244 flux ratio=4.384850649675 max_thermal=0.008465008044
+Iter number:7 k_eff=0.906834822807 flux ratio=4.384850652188 max_thermal=0.008465008093
+Iter number:8 k_eff=0.906834823192 flux ratio=4.384850654060 max_thermal=0.008465008128
+Iter number:9 k_eff=0.906834823457 flux ratio=4.384850655427 max_thermal=0.008465008153
+Iter number:10 k_eff=0.906834823641 flux ratio=4.384850656411 max_thermal=0.008465008170
+Iter number:11 k_eff=0.906834823769 flux ratio=4.384850657114 max_thermal=0.008465008182
+Iter number:12 k_eff=0.906834823858 flux ratio=4.384850657613 max_thermal=0.008465008191
+Iter number:13 k_eff=0.906834823919 flux ratio=4.384850657966 max_thermal=0.008465008197
+Iter number:14 k_eff=0.906834823962 flux ratio=4.384850658214 max_thermal=0.008465008201
+Iter number:15 k_eff=0.906834823992 flux ratio=4.384850658389 max_thermal=0.008465008204
+Iter number:16 k_eff=0.906834824013 flux ratio=4.384850658512 max_thermal=0.008465008206
+Iter number:17 k_eff=0.906834824028 flux ratio=4.384850658598 max_thermal=0.008465008207
+Iter number:18 k_eff=0.906834824038 flux ratio=4.384850658659 max_thermal=0.008465008208
+Iter number:19 k_eff=0.906834824046 flux ratio=4.384850658701 max_thermal=0.008465008209
+Iter number:20 k_eff=0.906834824051 flux ratio=4.384850658731 max_thermal=0.008465008209
+Iter number:21 k_eff=0.906834824054 flux ratio=4.384850658752 max_thermal=0.008465008210
+Iter number:22 k_eff=0.906834824057 flux ratio=4.384850658766 max_thermal=0.008465008210
+Iter number:23 k_eff=0.906834824058 flux ratio=4.384850658776 max_thermal=0.008465008210
+Iter number:24 k_eff=0.906834824059 flux ratio=4.384850658783 max_thermal=0.008465008210
+Iter number:25 k_eff=0.906834824060 flux ratio=4.384850658788 max_thermal=0.008465008210
+Cycle 7:
+   Numbers of active cells:       5434 21370 
+   Numbers of degrees of freedom: 23979 93275 
+Iter number:1 k_eff=0.906834695333 flux ratio=4.384848325238 max_thermal=0.008465012709
+Iter number:2 k_eff=0.906834781972 flux ratio=4.384848347736 max_thermal=0.008465013711
+Iter number:3 k_eff=0.906834782962 flux ratio=4.384848354071 max_thermal=0.008465013889
+Iter number:4 k_eff=0.906834784295 flux ratio=4.384848360506 max_thermal=0.008465014025
+Iter number:5 k_eff=0.906834785267 flux ratio=4.384848365788 max_thermal=0.008465014122
+Iter number:6 k_eff=0.906834785953 flux ratio=4.384848369732 max_thermal=0.008465014190
+Iter number:7 k_eff=0.906834786435 flux ratio=4.384848372568 max_thermal=0.008465014237
+Iter number:8 k_eff=0.906834786773 flux ratio=4.384848374574 max_thermal=0.008465014271
+Iter number:9 k_eff=0.906834787009 flux ratio=4.384848375984 max_thermal=0.008465014294
+Iter number:10 k_eff=0.906834787174 flux ratio=4.384848376972 max_thermal=0.008465014310
+Iter number:11 k_eff=0.906834787289 flux ratio=4.384848377663 max_thermal=0.008465014322
+Iter number:12 k_eff=0.906834787370 flux ratio=4.384848378146 max_thermal=0.008465014330
+Iter number:13 k_eff=0.906834787426 flux ratio=4.384848378484 max_thermal=0.008465014335
+Iter number:14 k_eff=0.906834787466 flux ratio=4.384848378720 max_thermal=0.008465014339
+Iter number:15 k_eff=0.906834787493 flux ratio=4.384848378885 max_thermal=0.008465014342
+Iter number:16 k_eff=0.906834787512 flux ratio=4.384848379000 max_thermal=0.008465014344
+Iter number:17 k_eff=0.906834787526 flux ratio=4.384848379081 max_thermal=0.008465014345
+Iter number:18 k_eff=0.906834787535 flux ratio=4.384848379137 max_thermal=0.008465014346
+Iter number:19 k_eff=0.906834787542 flux ratio=4.384848379176 max_thermal=0.008465014347
+Iter number:20 k_eff=0.906834787547 flux ratio=4.384848379204 max_thermal=0.008465014347
+Iter number:21 k_eff=0.906834787550 flux ratio=4.384848379223 max_thermal=0.008465014348
+Iter number:22 k_eff=0.906834787552 flux ratio=4.384848379237 max_thermal=0.008465014348
+Iter number:23 k_eff=0.906834787554 flux ratio=4.384848379246 max_thermal=0.008465014348
+Iter number:24 k_eff=0.906834787555 flux ratio=4.384848379253 max_thermal=0.008465014348
+Iter number:25 k_eff=0.906834787556 flux ratio=4.384848379257 max_thermal=0.008465014348
+Cycle 8:
+   Numbers of active cells:       6856 27001 
+   Numbers of degrees of freedom: 30285 117017 
+Iter number:1 k_eff=0.906834743244 flux ratio=4.384846479257 max_thermal=0.008465017253
+Iter number:2 k_eff=0.906834753823 flux ratio=4.384846519298 max_thermal=0.008465017557
+Iter number:3 k_eff=0.906834756845 flux ratio=4.384846524675 max_thermal=0.008465017731
+Iter number:4 k_eff=0.906834758500 flux ratio=4.384846529782 max_thermal=0.008465017851
+Iter number:5 k_eff=0.906834759516 flux ratio=4.384846533887 max_thermal=0.008465017934
+Iter number:6 k_eff=0.906834760177 flux ratio=4.384846537011 max_thermal=0.008465017993
+Iter number:7 k_eff=0.906834760619 flux ratio=4.384846539312 max_thermal=0.008465018034
+Iter number:8 k_eff=0.906834760921 flux ratio=4.384846540972 max_thermal=0.008465018062
+Iter number:9 k_eff=0.906834761129 flux ratio=4.384846542155 max_thermal=0.008465018083
+Iter number:10 k_eff=0.906834761274 flux ratio=4.384846542991 max_thermal=0.008465018097
+Iter number:11 k_eff=0.906834761374 flux ratio=4.384846543580 max_thermal=0.008465018106
+Iter number:12 k_eff=0.906834761444 flux ratio=4.384846543993 max_thermal=0.008465018113
+Iter number:13 k_eff=0.906834761493 flux ratio=4.384846544282 max_thermal=0.008465018118
+Iter number:14 k_eff=0.906834761527 flux ratio=4.384846544485 max_thermal=0.008465018121
+Iter number:15 k_eff=0.906834761551 flux ratio=4.384846544627 max_thermal=0.008465018124
+Iter number:16 k_eff=0.906834761567 flux ratio=4.384846544726 max_thermal=0.008465018125
+Iter number:17 k_eff=0.906834761579 flux ratio=4.384846544795 max_thermal=0.008465018127
+Iter number:18 k_eff=0.906834761587 flux ratio=4.384846544844 max_thermal=0.008465018127
+Iter number:19 k_eff=0.906834761593 flux ratio=4.384846544878 max_thermal=0.008465018128
+Iter number:20 k_eff=0.906834761597 flux ratio=4.384846544901 max_thermal=0.008465018128
+Iter number:21 k_eff=0.906834761599 flux ratio=4.384846544918 max_thermal=0.008465018129
+Iter number:22 k_eff=0.906834761601 flux ratio=4.384846544930 max_thermal=0.008465018129
+Iter number:23 k_eff=0.906834761603 flux ratio=4.384846544938 max_thermal=0.008465018129
+Iter number:24 k_eff=0.906834761604 flux ratio=4.384846544944 max_thermal=0.008465018129
+Iter number:25 k_eff=0.906834761604 flux ratio=4.384846544947 max_thermal=0.008465018129
+Cycle 9:
+   Numbers of active cells:       9166 35416 
+   Numbers of degrees of freedom: 40087 154355 
+Iter number:1 k_eff=0.906834746805 flux ratio=4.384846149270 max_thermal=0.008465018800
+Iter number:2 k_eff=0.906834739570 flux ratio=4.384846071221 max_thermal=0.008465018869
+Iter number:3 k_eff=0.906834742439 flux ratio=4.384846070361 max_thermal=0.008465018998
+Iter number:4 k_eff=0.906834743795 flux ratio=4.384846072627 max_thermal=0.008465019087
+Iter number:5 k_eff=0.906834744601 flux ratio=4.384846075162 max_thermal=0.008465019149
+Iter number:6 k_eff=0.906834745118 flux ratio=4.384846077315 max_thermal=0.008465019194
+Iter number:7 k_eff=0.906834745462 flux ratio=4.384846078976 max_thermal=0.008465019225
+Iter number:8 k_eff=0.906834745695 flux ratio=4.384846080201 max_thermal=0.008465019246
+Iter number:9 k_eff=0.906834745854 flux ratio=4.384846081084 max_thermal=0.008465019262
+Iter number:10 k_eff=0.906834745964 flux ratio=4.384846081712 max_thermal=0.008465019272
+Iter number:11 k_eff=0.906834746041 flux ratio=4.384846082157 max_thermal=0.008465019280
+Iter number:12 k_eff=0.906834746094 flux ratio=4.384846082469 max_thermal=0.008465019285
+Iter number:13 k_eff=0.906834746131 flux ratio=4.384846082689 max_thermal=0.008465019289
+Iter number:14 k_eff=0.906834746157 flux ratio=4.384846082842 max_thermal=0.008465019291
+Iter number:15 k_eff=0.906834746175 flux ratio=4.384846082950 max_thermal=0.008465019293
+Iter number:16 k_eff=0.906834746188 flux ratio=4.384846083025 max_thermal=0.008465019294
+Iter number:17 k_eff=0.906834746197 flux ratio=4.384846083078 max_thermal=0.008465019295
+Iter number:18 k_eff=0.906834746203 flux ratio=4.384846083115 max_thermal=0.008465019296
+Iter number:19 k_eff=0.906834746207 flux ratio=4.384846083141 max_thermal=0.008465019296
+Iter number:20 k_eff=0.906834746210 flux ratio=4.384846083159 max_thermal=0.008465019296
+Iter number:21 k_eff=0.906834746213 flux ratio=4.384846083171 max_thermal=0.008465019297
+Iter number:22 k_eff=0.906834746214 flux ratio=4.384846083180 max_thermal=0.008465019297
+Iter number:23 k_eff=0.906834746215 flux ratio=4.384846083186 max_thermal=0.008465019297
+Iter number:24 k_eff=0.906834746216 flux ratio=4.384846083191 max_thermal=0.008465019297
+Cycle 10:
+   Numbers of active cells:       10594 41197 
+   Numbers of degrees of freedom: 45467 179469 
+Iter number:1 k_eff=0.906780038935 flux ratio=4.384864081404 max_thermal=0.008464912703
+Iter number:2 k_eff=0.906833758951 flux ratio=4.384869346165 max_thermal=0.008465298273
+Iter number:3 k_eff=0.906834914205 flux ratio=4.384864425562 max_thermal=0.008465215350
+Iter number:4 k_eff=0.906835274876 flux ratio=4.384858287964 max_thermal=0.008465154357
+Iter number:5 k_eff=0.906835297124 flux ratio=4.384853804292 max_thermal=0.008465112059
+Iter number:6 k_eff=0.906835209298 flux ratio=4.384850953669 max_thermal=0.008465083180
+Iter number:7 k_eff=0.906835102972 flux ratio=4.384849195436 max_thermal=0.008465063473
+Iter number:8 k_eff=0.906835009160 flux ratio=4.384848100961 max_thermal=0.008465049973
+Iter number:9 k_eff=0.906834935068 flux ratio=4.384847404790 max_thermal=0.008465040684
+Iter number:10 k_eff=0.906834879514 flux ratio=4.384846951907 max_thermal=0.008465034267
+Iter number:11 k_eff=0.906834838999 flux ratio=4.384846651620 max_thermal=0.008465029821
+Iter number:12 k_eff=0.906834809916 flux ratio=4.384846449558 max_thermal=0.008465026733
+Iter number:13 k_eff=0.906834789237 flux ratio=4.384846312108 max_thermal=0.008465024585
+Iter number:14 k_eff=0.906834774618 flux ratio=4.384846217876 max_thermal=0.008465023089
+Iter number:15 k_eff=0.906834764321 flux ratio=4.384846152915 max_thermal=0.008465022045
+Iter number:16 k_eff=0.906834757086 flux ratio=4.384846107956 max_thermal=0.008465021317
+Iter number:17 k_eff=0.906834752010 flux ratio=4.384846076754 max_thermal=0.008465020809
+Iter number:18 k_eff=0.906834748452 flux ratio=4.384846055055 max_thermal=0.008465020454
+Iter number:19 k_eff=0.906834745960 flux ratio=4.384846039944 max_thermal=0.008465020207
+Iter number:20 k_eff=0.906834744216 flux ratio=4.384846029409 max_thermal=0.008465020033
+Iter number:21 k_eff=0.906834742995 flux ratio=4.384846022060 max_thermal=0.008465019912
+Iter number:22 k_eff=0.906834742141 flux ratio=4.384846016929 max_thermal=0.008465019828
+Iter number:23 k_eff=0.906834741544 flux ratio=4.384846013347 max_thermal=0.008465019768
+Iter number:24 k_eff=0.906834741126 flux ratio=4.384846010844 max_thermal=0.008465019727
+Iter number:25 k_eff=0.906834740834 flux ratio=4.384846009095 max_thermal=0.008465019698
+Iter number:26 k_eff=0.906834740629 flux ratio=4.384846007873 max_thermal=0.008465019678
+Iter number:27 k_eff=0.906834740486 flux ratio=4.384846007019 max_thermal=0.008465019664
+Iter number:28 k_eff=0.906834740386 flux ratio=4.384846006422 max_thermal=0.008465019654
+Iter number:29 k_eff=0.906834740316 flux ratio=4.384846006004 max_thermal=0.008465019647
+Iter number:30 k_eff=0.906834740267 flux ratio=4.384846005712 max_thermal=0.008465019642
+Iter number:31 k_eff=0.906834740233 flux ratio=4.384846005508 max_thermal=0.008465019639
+Iter number:32 k_eff=0.906834740209 flux ratio=4.384846005366 max_thermal=0.008465019637
+Iter number:33 k_eff=0.906834740193 flux ratio=4.384846005266 max_thermal=0.008465019635
+Iter number:34 k_eff=0.906834740181 flux ratio=4.384846005196 max_thermal=0.008465019634
+Iter number:35 k_eff=0.906834740173 flux ratio=4.384846005148 max_thermal=0.008465019633
+Iter number:36 k_eff=0.906834740167 flux ratio=4.384846005113 max_thermal=0.008465019632
+Iter number:37 k_eff=0.906834740163 flux ratio=4.384846005090 max_thermal=0.008465019632
+Iter number:38 k_eff=0.906834740160 flux ratio=4.384846005073 max_thermal=0.008465019632
+Iter number:39 k_eff=0.906834740158 flux ratio=4.384846005061 max_thermal=0.008465019632
+Iter number:40 k_eff=0.906834740157 flux ratio=4.384846005053 max_thermal=0.008465019631
+Iter number:41 k_eff=0.906834740156 flux ratio=4.384846005047 max_thermal=0.008465019631
+Iter number:42 k_eff=0.906834740155 flux ratio=4.384846005044 max_thermal=0.008465019631
+Cycle 11:
+   Numbers of active cells:       11749 47074 
+   Numbers of degrees of freedom: 50261 204523 
+Iter number:1 k_eff=0.906805395149 flux ratio=4.384872231023 max_thermal=0.008464861813
+Iter number:2 k_eff=0.906833353627 flux ratio=4.384863725577 max_thermal=0.008465049652
+Iter number:3 k_eff=0.906834008279 flux ratio=4.384854144970 max_thermal=0.008465027770
+Iter number:4 k_eff=0.906834374542 flux ratio=4.384849413955 max_thermal=0.008465020786
+Iter number:5 k_eff=0.906834559132 flux ratio=4.384847412067 max_thermal=0.008465018746
+Iter number:6 k_eff=0.906834648393 flux ratio=4.384846595090 max_thermal=0.008465018343
+Iter number:7 k_eff=0.906834691266 flux ratio=4.384846266355 max_thermal=0.008465018450
+Iter number:8 k_eff=0.906834712193 flux ratio=4.384846137013 max_thermal=0.008465018681
+Iter number:9 k_eff=0.906834722750 flux ratio=4.384846088708 max_thermal=0.008465018907
+Iter number:10 k_eff=0.906834728325 flux ratio=4.384846072840 max_thermal=0.008465019094
+Iter number:11 k_eff=0.906834731426 flux ratio=4.384846069459 max_thermal=0.008465019237
+Iter number:12 k_eff=0.906834733245 flux ratio=4.384846070461 max_thermal=0.008465019343
+Iter number:13 k_eff=0.906834734362 flux ratio=4.384846072648 max_thermal=0.008465019420
+Iter number:14 k_eff=0.906834735076 flux ratio=4.384846074842 max_thermal=0.008465019475
+Iter number:15 k_eff=0.906834735545 flux ratio=4.384846076677 max_thermal=0.008465019514
+Iter number:16 k_eff=0.906834735859 flux ratio=4.384846078098 max_thermal=0.008465019542
+Iter number:17 k_eff=0.906834736073 flux ratio=4.384846079156 max_thermal=0.008465019561
+Iter number:18 k_eff=0.906834736220 flux ratio=4.384846079925 max_thermal=0.008465019575
+Iter number:19 k_eff=0.906834736321 flux ratio=4.384846080478 max_thermal=0.008465019585
+Iter number:20 k_eff=0.906834736391 flux ratio=4.384846080872 max_thermal=0.008465019591
+Iter number:21 k_eff=0.906834736440 flux ratio=4.384846081150 max_thermal=0.008465019596
+Iter number:22 k_eff=0.906834736474 flux ratio=4.384846081347 max_thermal=0.008465019599
+Iter number:23 k_eff=0.906834736498 flux ratio=4.384846081485 max_thermal=0.008465019602
+Iter number:24 k_eff=0.906834736515 flux ratio=4.384846081582 max_thermal=0.008465019603
+Iter number:25 k_eff=0.906834736526 flux ratio=4.384846081650 max_thermal=0.008465019604
+Iter number:26 k_eff=0.906834736534 flux ratio=4.384846081698 max_thermal=0.008465019605
+Iter number:27 k_eff=0.906834736540 flux ratio=4.384846081731 max_thermal=0.008465019606
+Iter number:28 k_eff=0.906834736544 flux ratio=4.384846081755 max_thermal=0.008465019606
+Iter number:29 k_eff=0.906834736547 flux ratio=4.384846081771 max_thermal=0.008465019606
+Iter number:30 k_eff=0.906834736549 flux ratio=4.384846081783 max_thermal=0.008465019607
+Iter number:31 k_eff=0.906834736550 flux ratio=4.384846081791 max_thermal=0.008465019607
+Iter number:32 k_eff=0.906834736551 flux ratio=4.384846081796 max_thermal=0.008465019607
+Iter number:33 k_eff=0.906834736552 flux ratio=4.384846081800 max_thermal=0.008465019607
+@endcode
+
+We see that power iteration does converge faster after cycle 0 due to the initialization 
+with solution from last mesh iteration. 
+The contents of ``convergence_table'' are,
+@code
+0 4761 4761 0.906841960371 4.384056022582
+1 4761 10667 0.906837901031 4.385489087760
+2 4761 18805 0.906836075928 4.385466647499
+3 6629 27301 0.906835500111 4.385404580865
+4 12263 48095 0.906835001796 4.385381798734
+5 17501 69297 0.906834858174 4.384853823414
+6 19933 78605 0.906834824060 4.384850658788
+7 23979 93275 0.906834787556 4.384848379257
+8 30285 117017 0.906834761604 4.384846544947
+9 40087 154355 0.906834746216 4.384846083191
+10 45467 179469 0.906834740155 4.384846005044
+11 50261 204523 0.906834736552 4.384846081800
+@endcode
+The meanings of coloms are: number of mesh iteration, numbers of degrees of
+ freedom of fast energy group, numbers of DoFs of thermal group, converged 
+k-effective and the ratio between maximum of fast flux and maximum of thermal one.
+
+The grids of fast and thermal energy groups at mesh iteration #9 are shown 
+in following figure.
+
+@image html step-28.grid-0.9.order2.png
+@image html step-28.grid-1.9.order2.png
+
+We see that the grid of thermal group is much finner than the one of fast group. 
+The solutions on these grids are, (Note: flux are normalized with total fission
+source equal to 1)
+
+@image html step-28.solution-0.9.order2.png
+@image html step-28.solution-1.9.order2.png
+
+Then we plot the convergence data with polynomial order being equal to 1,2 and 3.
+
+@image html step-28.convergence.png
+
+The estimated ``exact'' k-effective = 0.906834721253 which is simply from last
+mesh iteration of polynomial order 3 minus 2e-10. We see that h-adaptive calculations
+deliver an algebraic convergence. And the higher polynomial order is, the faster mesh 
+iteration converges. In our problem, we need smaller number of DoFs to achieve same
+accuracy with higher polynoimal order.
diff --git a/deal.II/examples/step-28/doc/step-28.convergence.png b/deal.II/examples/step-28/doc/step-28.convergence.png
new file mode 100755 (executable)
index 0000000..65619bf
Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.convergence.png differ
diff --git a/deal.II/examples/step-28/doc/step-28.error-vs-dofs.png b/deal.II/examples/step-28/doc/step-28.error-vs-dofs.png
new file mode 100644 (file)
index 0000000..0240f37
Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.error-vs-dofs.png differ
diff --git a/deal.II/examples/step-28/doc/step-28.error-vs-time.png b/deal.II/examples/step-28/doc/step-28.error-vs-time.png
new file mode 100644 (file)
index 0000000..f3f8924
Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.error-vs-time.png differ
diff --git a/deal.II/examples/step-28/doc/step-28.grid-0.9.order2.png b/deal.II/examples/step-28/doc/step-28.grid-0.9.order2.png
new file mode 100755 (executable)
index 0000000..439a916
Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.grid-0.9.order2.png differ
diff --git a/deal.II/examples/step-28/doc/step-28.grid-1.9.order2.png b/deal.II/examples/step-28/doc/step-28.grid-1.9.order2.png
new file mode 100755 (executable)
index 0000000..7212932
Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.grid-1.9.order2.png differ
diff --git a/deal.II/examples/step-28/doc/step-28.solution-0.9.order2.png b/deal.II/examples/step-28/doc/step-28.solution-0.9.order2.png
new file mode 100755 (executable)
index 0000000..a01e238
Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.solution-0.9.order2.png differ
diff --git a/deal.II/examples/step-28/doc/step-28.solution-1.9.order2.png b/deal.II/examples/step-28/doc/step-28.solution-1.9.order2.png
new file mode 100755 (executable)
index 0000000..aa54318
Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.solution-1.9.order2.png differ
diff --git a/deal.II/examples/step-28/step-28.cc b/deal.II/examples/step-28/step-28.cc
new file mode 100644 (file)
index 0000000..304abc2
--- /dev/null
@@ -0,0 +1,2347 @@
+//TODO: remove direct (non-eigenvalue) problem possibility (remove
+//         ExtraneousSource and some stuff in run(), see isour)
+
+/*    $Id: step-6.cc,v 1.43 2006/02/09 02:54:11 wolf Exp $       */
+/*    Version: $Name:  $                                          */
+/*                                                                */
+/*    Copyright (C) 2006 by the deal.II authors and Yaqi Wang     */
+/*                                                                */
+/*    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.                        */
+
+                                 // @sect3{Include files}
+
+                                 // We start with a bunch of include
+                                 // files that have already been
+                                 // explained in previous tutorial
+                                 // programs:
+#include <base/timer.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/thread_management.h>
+#include <base/parameter_handler.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_refinement.h>
+#include <grid/grid_out.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+#include <iostream>
+
+#include <base/utilities.h>
+
+                                 // We use the next include file to
+                                 // access block vectors which provide
+                                 // us a convenient way to manage
+                                 // solution and right hand side
+                                 // vectors of all energy groups:
+#include <lac/block_vector.h>
+
+                                 // This include file is for
+                                 // transferring solutions from one
+                                 // mesh to another different mesh. We
+                                 // use it when we are initializing
+                                 // solutions after each mesh
+                                 // iteration:
+#include <numerics/solution_transfer.h>
+
+                                 // When integrating functions defined
+                                 // on one mesh against shape
+                                 // functions defined on a different
+                                 // mesh, we need a function @p
+                                 // get_finest_common_cells (as
+                                 // discussed in the introduction)
+                                 // which is defined in the following
+                                 // header file:
+#include <grid/grid_tools.h>
+
+                                 // Here are two more C++ standard
+                                 // headers that we use to define list
+                                 // data types as well as to fine-tune
+                                 // the output we generate:
+#include <list>
+#include <iomanip>
+
+                                // The last step is as in all
+                                // previous programs:
+using namespace dealii;
+
+
+                                // @sect3{Material data}
+
+                                 // First up, we need to define a
+                                 // class that provides material data
+                                 // (including diffusion coefficients,
+                                 // removal cross sections, scattering
+                                 // cross sections, fission cross
+                                 // sections and fission spectra) to
+                                 // the main class.
+                                 //
+                                 // The parameter to the constructor
+                                 // determines for how many energy
+                                 // groups we set up the relevant
+                                 // tables. At present, this program
+                                 // only includes data for 2 energy
+                                 // groups, but a more sophisticated
+                                 // program may be able to initialize
+                                 // the data structures for more
+                                 // groups as well, depending on how
+                                 // many energy groups are selected in
+                                 // the parameter file.
+                                 //
+                                 // For each of the different
+                                 // coefficient types, there is one
+                                 // function that returns the value of
+                                 // this coefficient for a particular
+                                 // energy group (or combination of
+                                 // energy groups, as for the
+                                 // distribution cross section
+                                 // $\chi_g\nu\Sigma_{f,g'}$ or
+                                 // scattering cross section
+                                 // $\Sigma_{s,g'\to g}$). In addition
+                                 // to the energy group or groups,
+                                 // these coefficients depend on the
+                                 // type of fuel or control rod, as
+                                 // explained in the introduction. The
+                                 // functions therefore take an
+                                 // additional parameter, @p
+                                 // material_id, that identifies the
+                                 // particular kind of rod. Within
+                                 // this program, we use
+                                 // <code>n_materials=8</code>
+                                 // different kinds of rods.
+                                 //
+                                 // Except for the scattering cross
+                                 // section, each of the coefficients
+                                 // therefore can be represented as an
+                                 // entry in a two-dimensional array
+                                 // of floating point values indexed
+                                 // by the energy group number as well
+                                 // as the material ID. The Table
+                                 // class template is the ideal way to
+                                 // store such data. Finally, the
+                                 // scattering coefficient depends on
+                                 // both two energy group indices and
+                                 // therefore needs to be stored in a
+                                 // three-dimensional array, for which
+                                 // we again use the Table class,
+                                 // where this time the first template
+                                 // argument (denoting the
+                                 // dimensionality of the array) of
+                                 // course needs to be three:
+class MaterialData
+{
+  public:
+    MaterialData (const unsigned int n_groups);
+    double get_diffusion_coefficient (const unsigned int group,
+                                     const unsigned int material_id) const;
+    double get_removal_XS (const unsigned int group,
+                          const unsigned int material_id) const;
+    double get_fission_XS (const unsigned int group,
+                          const unsigned int material_id) const;
+    double get_fission_dist_XS (const unsigned int group_1,
+                               const unsigned int group_2,
+                               const unsigned int material_id) const;
+    double get_scattering_XS (const unsigned int group_1,
+                             const unsigned int group_2,
+                             const unsigned int material_id) const;
+    double get_fission_spectrum (const unsigned int group,
+                                const unsigned int material_id) const;
+    
+  private:
+    const unsigned int n_groups;
+    const unsigned int n_materials;
+
+    Table<2,double> diffusion;
+    Table<2,double> sigma_r;
+    Table<2,double> nu_sigma_f;
+    Table<3,double> sigma_s;
+    Table<2,double> chi;
+};
+
+                                 // The constructor of the class is
+                                 // used to initialize all the
+                                 // material data arrays. It takes the
+                                 // number of energy groups as an
+                                 // argument (an throws an error if
+                                 // that value is not equal to two,
+                                 // since at presently only data for
+                                 // two energy groups is implemented;
+                                 // however, using this, the function
+                                 // remains flexible and extendible
+                                 // into the future). In the member
+                                 // initialization part at the
+                                 // beginning, it also resizes the
+                                 // arrays to their correct sizes.
+                                 //
+                                 // At present, material data is
+                                 // stored for 8 different types of
+                                 // material. This, as well, may
+                                 // easily be extended in the future.
+MaterialData::MaterialData (const unsigned int n_groups)
+               :
+               n_groups (n_groups),
+               n_materials (8),
+               diffusion (n_materials, n_groups),
+               sigma_r (n_materials, n_groups),
+               nu_sigma_f (n_materials, n_groups),
+               sigma_s (n_materials, n_groups, n_groups),
+               chi (n_materials, n_groups)
+{
+  switch (n_groups)
+    {
+      case 2:
+      {
+        for (unsigned int m=0; m<n_materials; ++m)
+          {
+            diffusion[m][0] = 1.2;
+            diffusion[m][1] = 0.4;
+            chi[m][0]       = 1.0;
+            chi[m][1]       = 0.0;
+            sigma_r[m][0]   = 0.03;
+            for (unsigned int group_1=0; group_1<n_groups; ++group_1)
+              for (unsigned int group_2=0; group_2<n_groups; ++ group_2)
+                sigma_s[m][group_1][group_2]   = 0.0;
+          }
+    
+  
+        diffusion[5][1]  = 0.2;
+  
+        sigma_r[4][0]    = 0.026;
+        sigma_r[5][0]    = 0.051;
+        sigma_r[6][0]    = 0.026;
+        sigma_r[7][0]    = 0.050;
+  
+        sigma_r[0][1]    = 0.100;
+        sigma_r[1][1]    = 0.200;
+        sigma_r[2][1]    = 0.250;
+        sigma_r[3][1]    = 0.300;
+        sigma_r[4][1]    = 0.020;
+        sigma_r[5][1]    = 0.040;
+        sigma_r[6][1]    = 0.020;
+        sigma_r[7][1]    = 0.800;
+  
+        nu_sigma_f[0][0] = 0.0050;
+        nu_sigma_f[1][0] = 0.0075;
+        nu_sigma_f[2][0] = 0.0075;
+        nu_sigma_f[3][0] = 0.0075;
+        nu_sigma_f[4][0] = 0.000;
+        nu_sigma_f[5][0] = 0.000;
+        nu_sigma_f[6][0] = 1e-7;
+        nu_sigma_f[7][0] = 0.00;
+  
+        nu_sigma_f[0][1] = 0.125;
+        nu_sigma_f[1][1] = 0.300;
+        nu_sigma_f[2][1] = 0.375;
+        nu_sigma_f[3][1] = 0.450;
+        nu_sigma_f[4][1] = 0.000;
+        nu_sigma_f[5][1] = 0.000;
+        nu_sigma_f[6][1] = 3e-6;
+        nu_sigma_f[7][1] = 0.00;
+    
+        sigma_s[0][0][1] = 0.020;
+        sigma_s[1][0][1] = 0.015;
+        sigma_s[2][0][1] = 0.015;
+        sigma_s[3][0][1] = 0.015;
+        sigma_s[4][0][1] = 0.025;
+        sigma_s[5][0][1] = 0.050;
+        sigma_s[6][0][1] = 0.025;
+        sigma_s[7][0][1] = 0.010;
+
+        break;
+      }
+
+
+      default:
+            Assert (false,
+                    ExcMessage ("Presently, only data for 2 groups is implemented"));
+    }          
+}
+
+
+                                 // Next are the functions that return
+                                 // the coefficient values for given
+                                 // materials and energy groups. All
+                                 // they do is to make sure that the
+                                 // given arguments are within the
+                                 // allowed ranges, and then look the
+                                 // respective value up in the
+                                 // corresponding tables:
+double
+MaterialData::get_diffusion_coefficient (const unsigned int group,
+                                        const unsigned int material_id) const
+{
+  Assert (group < n_groups, 
+         ExcIndexRange (group, 0, n_groups));
+  Assert (material_id < n_materials, 
+         ExcIndexRange (material_id, 0, n_materials));
+
+  return diffusion[material_id][group];
+}
+
+
+
+double
+MaterialData::get_removal_XS (const unsigned int group,
+                              const unsigned int material_id) const
+{
+  Assert (group < n_groups, 
+         ExcIndexRange (group, 0, n_groups));
+  Assert (material_id < n_materials, 
+         ExcIndexRange (material_id, 0, n_materials));
+
+  return sigma_r[material_id][group];
+}
+
+
+double
+MaterialData::get_fission_XS (const unsigned int group,
+                              const unsigned int material_id) const
+{
+  Assert (group < n_groups, 
+         ExcIndexRange (group, 0, n_groups));
+  Assert (material_id < n_materials, 
+         ExcIndexRange (material_id, 0, n_materials));
+
+  return nu_sigma_f[material_id][group];
+}
+
+
+
+double
+MaterialData::get_scattering_XS (const unsigned int group_1,
+                                 const unsigned int group_2,
+                                 const unsigned int material_id) const
+{
+  Assert (group_1 < n_groups, 
+         ExcIndexRange (group_1, 0, n_groups));
+  Assert (group_2 < n_groups, 
+         ExcIndexRange (group_2, 0, n_groups));
+  Assert (material_id < n_materials, 
+         ExcIndexRange (material_id, 0, n_materials));
+
+  return sigma_s[material_id][group_1][group_2];
+}
+
+
+
+double
+MaterialData::get_fission_spectrum (const unsigned int group,
+                                    const unsigned int material_id) const 
+{
+  Assert (group < n_groups, 
+         ExcIndexRange (group, 0, n_groups));
+  Assert (material_id < n_materials, 
+         ExcIndexRange (material_id, 0, n_materials));
+
+  return chi[material_id][group];
+}
+
+
+                                 // The function computing the fission
+                                 // distribution cross section is
+                                 // slightly different, since it
+                                 // computes its value as the product
+                                 // of two other coefficients. We
+                                 // don't need to check arguments
+                                 // here, since this already happens
+                                 // when we call the two other
+                                 // functions involved, even though it
+                                 // would probably not hurt either:
+double
+MaterialData::get_fission_dist_XS (const unsigned int group_1,
+                                   const unsigned int group_2,
+                                   const unsigned int material_id) const 
+{
+  return (get_fission_spectrum(group_1, material_id) *
+         get_fission_XS(group_2, material_id));
+}
+
+
+
+                                // @sect3{The <code>EnergyGroup</code> class}
+
+                                // The first interesting class is the
+                                // one that contains everything that
+                                // is specific to a single energy
+                                // group. To group things that belong
+                                // together into individual objects,
+                                // we declare a structure that holds
+                                // the Triangulation and DoFHandler
+                                // objects for the mesh used for a
+                                // single energy group, and a number
+                                // of other objects and member
+                                // functions that we will discuss in
+                                // the following sections.
+                                //
+                                // The main reason for this class is
+                                // as follows: for both the forward
+                                // problem (with a specified right
+                                // hand side) as well as for the
+                                // eigenvalue problem, one typically
+                                // solves a sequence of problems for
+                                // a single energy group each, rather
+                                // than the fully coupled
+                                // problem. This becomes
+                                // understandable once one realizes
+                                // that the system matrix for a
+                                // single energy group is symmetric
+                                // and positive definite (it is
+                                // simply a diffusion operator),
+                                // whereas the matrix for the fully
+                                // coupled problem is generally
+                                // nonsymmetric and not definite. It
+                                // is also very large and quite full
+                                // if more than a few energy groups
+                                // are involved.
+                                //
+                                // Let us first look at the equation
+                                // to solve in the case of an
+                                // external right hand side (for the time
+                                // independent case):
+                                // @f{eqnarray*}
+                                // -\nabla \cdot(D_g(x) \nabla \phi_g(x))
+                                // +
+                                // \Sigma_{r,g}(x)\phi_g(x)
+                                // =
+                                // \chi_g\sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}(x)
+                                // +
+                                // \sum_{g'\ne g}\Sigma_{s,g'\to g}(x)\phi_{g'}(x)
+                                // +
+                                // s_{\mathrm{ext},g}(x)
+                                // @f}
+                                //
+                                // We would typically solve this
+                                // equation by moving all the terms
+                                // on the right hand side with $g'=g$
+                                // to the left hand side, and solving
+                                // for $\phi_g$. Of course, we don't
+                                // know $\phi_{g'}$ yet, since the
+                                // equations for those variables
+                                // include right hand side terms
+                                // involving $\phi_g$. What one
+                                // typically does in such situations
+                                // is to iterate: compute 
+                                // @f{eqnarray*}
+                                // -\nabla \cdot(D_g(x) \nabla \phi^{(n)}_g(x))
+                                // &+&
+                                // \Sigma_{r,g}(x)\phi^{(n)}_g(x)
+                                // \\ &=&
+                                // \chi_g\sum_{g'=1}^{g-1}\nu\Sigma_{f,g'}(x)\phi^{(n)}_{g'}(x)
+                                // +
+                                // \chi_g\sum_{g'=g}^G\nu\Sigma_{f,g'}(x)\phi^{(n-1)}_{g'}(x)
+                                // +
+                                // \sum_{g'\ne g, g'<g}\Sigma_{s,g'\to g}(x)\phi^{(n)}_{g'}(x)
+                                // +
+                                // \sum_{g'\ne g, g'>g}\Sigma_{s,g'\to g}(x)\phi^{(n-1)}_{g'}(x)
+                                // +
+                                // s_{\mathrm{ext},g}(x)
+                                // @f}
+                                //
+                                // In other words, we solve the
+                                // equation one by one, using values
+                                // for $\phi_{g'}$ from the previous
+                                // iteration $n-1$ if $g'\ge g$ and
+                                // already computed values for
+                                // $\phi_{g'}$ from the present
+                                // iteration if $g'<g$.
+                                //
+                                // When computing the eigenvalue, we
+                                // do a very similar iteration,
+                                // except that we have no external
+                                // right hand side and that the
+                                // solution is scaled after each
+                                // iteration as explained in the
+                                // introduction.
+                                //
+                                // In either case, these two cases
+                                // can be treated jointly if all we
+                                // do is to equip the following class
+                                // with these abilities: (i) form the
+                                // left hand side matrix, (ii) form
+                                // the in-group right hand side
+                                // contribution, i.e. involving the
+                                // extraneous source, and (iii) form
+                                // that contribution to the right
+                                // hand side that stems from group
+                                // $g'$. This class does exactly
+                                // these tasks (as well as some
+                                // book-keeping, such as mesh
+                                // refinement, setting up matrices
+                                // and vectors, etc). On the other
+                                // hand, the class itself has no idea
+                                // how many energy groups there are,
+                                // and in particular how they
+                                // interact, i.e. the decision of how
+                                // the outer iteration looks (and
+                                // consequently whether we solve an
+                                // eigenvalue or a direct problem) is
+                                // left to the
+                                // NeutronDiffusionProblem class
+                                // further down below in this
+                                // program.
+                                //
+                                // So let us go through the class and
+                                // its interface:
+template <int dim>
+class EnergyGroup
+{
+  public:
+
+                                    // @sect5{Public member functions}
+                                    //
+                                    // The class has a good number of
+                                    // public member functions, since
+                                    // its the way it operates is
+                                    // controlled from the outside,
+                                    // and therefore all functions
+                                    // that do something significant
+                                    // need to be called from another
+                                    // class. Let's start off with
+                                    // book-keeping: the class
+                                    // obviously needs to know which
+                                    // energy group it represents,
+                                    // which material data to use,
+                                    // and from what coarse grid to
+                                    // start. The constructor takes
+                                    // this information and
+                                    // initializes the relevant
+                                    // member variables with that
+                                    // (see below).
+                                    //
+                                    // Then we also need functions
+                                    // that set up the linear system,
+                                    // i.e. correctly size the matrix
+                                    // and its sparsity pattern, etc,
+                                    // given a finite element object
+                                    // to use. The
+                                    // <code>setup_linear_system</code>
+                                    // function does that. Finally,
+                                    // for this initial block, there
+                                    // are two functions that return
+                                    // the number of active cells and
+                                    // degrees of freedom used in
+                                    // this object -- using this, we
+                                    // can make the triangulation and
+                                    // DoF handler member variables
+                                    // private, and do not have to
+                                    // grant external use to it,
+                                    // enhancing encapsulation:
+    EnergyGroup (const unsigned int        group,
+                const MaterialData       &material_data,
+                const Triangulation<dim> &coarse_grid,
+                 const FiniteElement<dim> &fe);
+
+    void setup_linear_system ();
+
+    unsigned int n_active_cells () const;
+    unsigned int n_dofs () const;
+    
+                                     // Then there are functions that
+                                     // assemble the linear system for
+                                     // each iteration and the present
+                                     // energy group. Note that the
+                                     // matrix is independent of the
+                                     // iteration number, so only has
+                                     // to be computed once for each
+                                     // refinement cycle. The
+                                     // situation is a bit more
+                                     // involved for the right hand
+                                     // side that has to be updated in
+                                     // each inverse power iteration,
+                                     // and that is further
+                                     // complicated by the fact that
+                                     // computing it may involve
+                                     // several different meshes as
+                                     // explained in the
+                                     // introduction. To make things
+                                     // more flexible with regard to
+                                     // solving the forward or the
+                                     // eigenvalue problem, we split
+                                     // the computation of the right
+                                     // hand side into a function that
+                                     // assembles the extraneous
+                                     // source and in-group
+                                     // contributions (which we will
+                                     // call with a zero function as
+                                     // source terms for the
+                                     // eigenvalue problem) and one
+                                     // that computes contributions to
+                                     // the right hand side from
+                                     // another energy group:
+    void assemble_system_matrix ();
+    void assemble_ingroup_rhs (const Function<dim> &extraneous_source);
+    void assemble_cross_group_rhs (const EnergyGroup<dim> &g_prime);
+
+                                    // Next we need a set of
+                                    // functions that actually
+                                    // compute the solution of a
+                                    // linear system, and do
+                                    // something with it (such as
+                                    // computing the fission source
+                                    // contribution mentioned in the
+                                    // introduction, writing
+                                    // graphical information to an
+                                    // output file, computing error
+                                    // indicators, or actually
+                                    // refining the grid based on
+                                    // these criteria and thresholds
+                                    // for refinement and
+                                    // coarsening). All these
+                                    // functions will later be called
+                                    // from the driver class
+                                    // <code>NeutronDiffusionProblem</code>,
+                                    // or any other class you may
+                                    // want to implement to solve a
+                                    // problem involving the neutron
+                                    // flux equations:
+    void   solve ();
+    
+    double get_fission_source () const;
+    
+    void   output_results (const unsigned int cycle) const;
+    
+    void   estimate_errors (Vector<float> &error_indicators) const;
+    
+    void   refine_grid (const Vector<float> &error_indicators,
+                        const double         refine_threshold,
+                        const double         coarsen_threshold);
+
+                                    // @sect5{Public data members}
+                                    //
+                                    // As is good practice in object
+                                    // oriented programming, we hide
+                                    // most data members by making
+                                    // them private. However, we have
+                                    // to grant the class that drives
+                                    // the process access to the
+                                    // solution vector as well as the
+                                    // solution of the previous
+                                    // iteration, since in the power
+                                    // iteration, the solution vector
+                                    // is scaled in every iteration
+                                    // by the present guess of the
+                                    // eigenvalue we are looking for:
+  public:
+    
+    Vector<double> solution;
+    Vector<double> solution_old;
+
+
+                                    // @sect5{Private data members}
+                                    //
+                                     // The rest of the data members
+                                     // are private. Compared to all
+                                     // the previous tutorial
+                                     // programs, the only new data
+                                     // members are an integer storing
+                                     // which energy group this object
+                                     // represents, and a reference to
+                                     // the material data object that
+                                     // this object's constructor gets
+                                     // passed from the driver
+                                     // class. Likewise, the
+                                     // constructor gets a reference
+                                     // to the finite element object
+                                     // we are to use.
+                                     //
+                                     // Finally, we have to apply
+                                     // boundary values to the linear
+                                     // system in each iteration,
+                                     // i.e. quite frequently. Rather
+                                     // than interpolating them every
+                                     // time, we interpolate them once
+                                     // on each new mesh and then
+                                     // store them along with all the
+                                     // other data of this class:
+  private:
+    
+    const unsigned int            group;
+    const MaterialData           &material_data;
+        
+    Triangulation<dim>            triangulation;
+    const FiniteElement<dim>     &fe;
+    DoFHandler<dim>               dof_handler;
+
+    SparsityPattern               sparsity_pattern;
+    SparseMatrix<double>          system_matrix;
+    
+    Vector<double>                system_rhs;
+    
+    std::map<unsigned int,double> boundary_values;
+    ConstraintMatrix              hanging_node_constraints;
+    
+
+                                    // @sect5{Private member functionss}
+                                    //
+                                     // There is one private member
+                                     // function in this class. It
+                                     // recursively walks over cells
+                                     // of two meshes to compute the
+                                     // cross-group right hand side
+                                     // terms. The algorithm for this
+                                     // is explained in the
+                                     // introduction to this
+                                     // program. The arguments to this
+                                     // function are a reference to an
+                                     // object representing the energy
+                                     // group against which we want to
+                                     // integrate a right hand side
+                                     // term, an iterator to a cell of
+                                     // the mesh used for the present
+                                     // energy group, an iterator to a
+                                     // corresponding cell on the
+                                     // other mesh, and the matrix
+                                     // that interpolates the degrees
+                                     // of freedom from the coarser of
+                                     // the two cells to the finer
+                                     // one:
+  private:
+
+    void
+    assemble_cross_group_rhs_recursive (const EnergyGroup<dim>                        &g_prime,
+                                        const typename DoFHandler<dim>::cell_iterator &cell_g,
+                                        const typename DoFHandler<dim>::cell_iterator &cell_g_prime,
+                                        const FullMatrix<double>                       prolongation_matrix);
+};
+
+
+                                 // @sect4{Implementation of the <code>EnergyGroup</code> class}
+
+                                 // The first few functions of this
+                                 // class are mostly
+                                 // self-explanatory. The constructor
+                                 // only sets a few data members and
+                                 // creates a copy of the given
+                                 // triangulation as the base for the
+                                 // triangulation used for this energy
+                                 // group. The next two functions
+                                 // simply return data from private
+                                 // data members, thereby enabling us
+                                 // to make these data members
+                                 // private.
+template <int dim>
+EnergyGroup<dim>::EnergyGroup (const unsigned int        group,
+                              const MaterialData       &material_data,
+                              const Triangulation<dim> &coarse_grid,
+                               const FiniteElement<dim> &fe)
+               :
+                group (group),
+               material_data (material_data),
+                fe (fe),
+               dof_handler (triangulation)
+{
+  triangulation.copy_triangulation (coarse_grid);
+  dof_handler.distribute_dofs (fe);
+}
+
+
+
+template <int dim>
+unsigned int
+EnergyGroup<dim>::n_active_cells () const
+{
+  return triangulation.n_active_cells ();
+}
+
+
+
+template <int dim>
+unsigned int
+EnergyGroup<dim>::n_dofs () const
+{
+  return dof_handler.n_dofs ();
+}
+
+
+
+                                 // @sect5{<code>EnergyGroup::setup_linear_system</code>}
+                                 //
+                                 // The first "real" function is the
+                                 // one that sets up the mesh,
+                                 // matrices, etc, on the new mesh or
+                                 // after mesh refinement. We use this
+                                 // function to initialize sparse
+                                 // system matrices, and the right
+                                 // hand side vector. If the solution
+                                 // vector has never been set before
+                                 // (as indicated by a zero size), we
+                                 // also initialize it and set it to a
+                                 // default value. We don't do that if
+                                 // it already has a non-zero size
+                                 // (i.e. this function is called
+                                 // after mesh refinement) since in
+                                 // that case we want to preserve the
+                                 // solution across mesh refinement
+                                 // (something we do in the
+                                 // <code>EnergyGroup::refine_grid</code>
+                                 // function).
+template <int dim>
+void
+EnergyGroup<dim>::setup_linear_system ()
+{
+  const unsigned int n_dofs = dof_handler.n_dofs();
+      
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                           hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  system_matrix.clear ();
+      
+  sparsity_pattern.reinit (n_dofs, n_dofs,
+                           dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  hanging_node_constraints.condense (sparsity_pattern);
+  sparsity_pattern.compress ();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  system_rhs.reinit (n_dofs);
+
+  if (solution.size() == 0)
+    {
+      solution.reinit (n_dofs);
+      solution_old.reinit(n_dofs);
+      solution_old = 1.0;
+      solution = solution_old;
+    }
+
+
+                                   // At the end of this function, we
+                                   // update the list of boundary
+                                   // nodes and their values, by first
+                                   // clearing this list and the
+                                   // re-interpolating boundary values
+                                   // (remember that this function is
+                                   // called after first setting up
+                                   // the mesh, and each time after
+                                   // mesh refinement).
+                                   //
+                                   // To understand the code, it is
+                                   // necessary to realize that we
+                                   // create the mesh using the
+                                   // <code>GridGenerator::subdivided_hyper_rectangle</code>
+                                   // function (in
+                                   // <code>NeutronDiffusionProblem::initialize_problem</code>)
+                                   // where we set the last parameter
+                                   // to <code>true</code>. This means that
+                                   // boundaries of the domain are
+                                   // "colored", i.e. the four (or
+                                   // six, in 3d) sides of the domain
+                                   // are assigned different boundary
+                                   // indicators. As it turns out, the
+                                   // bottom boundary gets indicator
+                                   // zero, the top one boundary
+                                   // indicator one, and left and
+                                   // right boundaries get indicators
+                                   // two and three, respectively.
+                                   //
+                                   // In this program, we simulate
+                                   // only one, namely the top right,
+                                   // quarter of a reactor. That is,
+                                   // we want to interpolate boundary
+                                   // conditions only on the top and
+                                   // right boundaries, while do
+                                   // nothing on the bottom and left
+                                   // boundaries (i.e. impose natural,
+                                   // no-flux Neumann boundary
+                                   // conditions). This is most easily
+                                   // generalized to arbitrary
+                                   // dimension by saying that we want
+                                   // to interpolate on those
+                                   // boundaries with indicators 1, 3,
+                                   // ..., which we do in the
+                                   // following loop (note that calls
+                                   // to
+                                   // <code>VectorTools::interpolate_boundary_values</code>
+                                   // are additive, i.e. they do not
+                                   // first clear the boundary value
+                                   // map):
+  boundary_values.clear();
+    
+  for (unsigned int i=0; i<dim; ++i)
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                              2*i+1,
+                                              ZeroFunction<dim>(),
+                                              boundary_values);
+}
+
+
+
+                                 // @sect5{<code>EnergyGroup::assemble_system_matrix</code>}
+                                 //
+                                 // Next we need functions assembling
+                                 // the system matrix and right hand
+                                 // sides. Assembling the matrix is
+                                 // straightforward given the
+                                 // equations outlined in the
+                                 // introduction as well as what we've
+                                 // seen in previous example
+                                 // programs. Note the use of
+                                 // <code>cell->material_id()</code> to get at
+                                 // the kind of material from which a
+                                 // cell is made up of. Note also how
+                                 // we set the order of the quadrature
+                                 // formula so that it is always
+                                 // appropriate for the finite element
+                                 // in use.
+                                 //
+                                 // Finally, note that since we only
+                                 // assemble the system matrix here,
+                                 // we can't yet eliminate boundary
+                                 // values (we need the right hand
+                                 // side vector for this). We defer
+                                 // this to the <code>EnergyGroup::solve</code>
+                                 // function, at which point all the
+                                 // information is available.
+template <int dim>
+void
+EnergyGroup<dim>::assemble_system_matrix () 
+{
+  const QGauss<dim>  quadrature_formula(fe.degree + 1);
+      
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values    |  update_gradients |
+                           update_JxW_values);
+      
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.n_quadrature_points;
+      
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs (dofs_per_cell);
+      
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+      
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+      
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+         
+      fe_values.reinit (cell);
+
+      const double diffusion_coefficient
+       = material_data.get_diffusion_coefficient (group, cell->material_id());
+      const double removal_XS
+       = material_data.get_removal_XS (group,cell->material_id());
+         
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+            cell_matrix(i,j) += ((diffusion_coefficient *
+                                  fe_values.shape_grad(i,q_point) *
+                                  fe_values.shape_grad(j,q_point)
+                                  +
+                                  removal_XS *
+                                  fe_values.shape_value(i,q_point) *
+                                  fe_values.shape_value(j,q_point))
+                                 *
+                                 fe_values.JxW(q_point));
+      
+      cell->get_dof_indices (local_dof_indices);
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         system_matrix.add (local_dof_indices[i],
+                            local_dof_indices[j],
+                            cell_matrix(i,j));
+    }
+  
+  hanging_node_constraints.condense (system_matrix);
+}
+
+
+
+                                 // @sect5{<code>EnergyGroup::assemble_ingroup_rhs</code>}
+                                 //
+                                 // As explained in the documentation
+                                 // of the <code>EnergyGroup</code> class, we
+                                 // split assembling the right hand
+                                 // side into two parts: the ingroup
+                                 // and the cross-group
+                                 // couplings. First, we need a
+                                 // function to assemble the right
+                                 // hand side of one specific group
+                                 // here, i.e. including an extraneous
+                                 // source (that we will set to zero
+                                 // for the eigenvalue problem) as
+                                 // well as the ingroup fission
+                                 // contributions.  (In-group
+                                 // scattering has already been
+                                 // accounted for with the definition
+                                 // of removal cross section.) The
+                                 // function's workings are pretty
+                                 // standard as far as assembling
+                                 // right hand sides go, and therefore
+                                 // does not require more comments
+                                 // except that we mention that the
+                                 // right hand side vector is set to
+                                 // zero at the beginning of the
+                                 // function -- something we are not
+                                 // going to do for the cross-group
+                                 // terms that simply add to the right
+                                 // hand side vector.
+template <int dim>
+void EnergyGroup<dim>::assemble_ingroup_rhs (const Function<dim> &extraneous_source)
+{
+  system_rhs.reinit (dof_handler.n_dofs());
+
+  const QGauss<dim>  quadrature_formula (fe.degree + 1);
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+  
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values    |  update_q_points  |
+                           update_JxW_values);
+  
+  Vector<double>            cell_rhs (dofs_per_cell);
+  std::vector<double>       extraneous_source_values (n_q_points);
+  std::vector<double>       solution_old_values (n_q_points);
+  
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  
+  for (; cell!=endc; ++cell)
+    {
+      cell_rhs = 0;
+      
+      fe_values.reinit (cell);
+      
+      const double fission_dist_XS
+       = material_data.get_fission_dist_XS (group, group, cell->material_id());
+
+      extraneous_source.value_list (fe_values.get_quadrature_points(),
+                                   extraneous_source_values);
+
+      fe_values.get_function_values (solution_old, solution_old_values);
+      
+      cell->get_dof_indices (local_dof_indices);
+      
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         cell_rhs(i) += ((extraneous_source_values[q_point]
+                           +
+                          fission_dist_XS *
+                           solution_old_values[q_point]) *
+                         fe_values.shape_value(i,q_point) *
+                         fe_values.JxW(q_point));
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       system_rhs(local_dof_indices[i]) += cell_rhs(i);      
+    }
+}
+
+
+
+                                 // @sect5{<code>EnergyGroup::assemble_cross_group_rhs</code>}
+                                 //
+                                 // The more interesting function for
+                                 // assembling the right hand side
+                                 // vector for the equation of a
+                                 // single energy group is the one
+                                 // that couples energy group $g$ and
+                                 // $g'$. As explained in the
+                                 // introduction, we first have to
+                                 // find the set of cells common to
+                                 // the meshes of the two energy
+                                 // groups. First we call
+                                 // <code>get_finest_common_cells</code> to
+                                 // obtain this list of pairs of
+                                 // common cells from both
+                                 // meshes. Both cells in a pair may
+                                 // not be active but at least one of
+                                 // them is. We then hand each of
+                                 // these cell pairs off to a function
+                                 // tha computes the right hand side
+                                 // terms recursively.
+                                 //
+                                 // Note that ingroup coupling is
+                                 // handled already before, so we exit
+                                 // the function early if $g=g'$.
+template <int dim>
+void EnergyGroup<dim>::assemble_cross_group_rhs (const EnergyGroup<dim> &g_prime)
+{
+  if (group == g_prime.group)
+    return;
+  
+  const std::list<std::pair<typename DoFHandler<dim>::cell_iterator, 
+    typename DoFHandler<dim>::cell_iterator> >
+    cell_list
+    = GridTools::get_finest_common_cells (dof_handler, 
+                                         g_prime.dof_handler);
+      
+  typename std::list<std::pair<typename DoFHandler<dim>::cell_iterator, 
+                               typename DoFHandler<dim>::cell_iterator> >
+    ::const_iterator 
+    cell_iter = cell_list.begin();
+      
+  for (; cell_iter!=cell_list.end(); ++cell_iter)
+    {
+      FullMatrix<double> unit_matrix (fe.dofs_per_cell);
+      for (unsigned int i=0; i<unit_matrix.m(); ++i)
+       unit_matrix(i,i) = 1;
+      assemble_cross_group_rhs_recursive (g_prime,
+                                          cell_iter->first,
+                                         cell_iter->second,
+                                         unit_matrix);
+    }
+}
+
+
+
+                                 // @sect5{<code>EnergyGroup::assemble_cross_group_rhs_recursive</code>}
+                                 //
+                                 // This is finally the function that
+                                 // handles assembling right hand side
+                                 // terms on potentially different
+                                 // meshes recursively, using the
+                                 // algorithm described in the
+                                 // introduction. The function takes a
+                                 // reference to the object
+                                 // representing energy group $g'$, as
+                                 // well as iterators to corresponding
+                                 // cells in the meshes for energy
+                                 // groups $g$ and $g'$. At first,
+                                 // i.e. when this function is called
+                                 // from the one above, these two
+                                 // cells will be matching cells on
+                                 // two meshes; however, one of the
+                                 // two may be further refined, and we
+                                 // will call the function recursively
+                                 // with one of the two iterators
+                                 // replaced by one of the children of
+                                 // the original cell.
+                                 //
+                                 // The last argument is the matrix
+                                 // product matrix $B_{c^{(k)}}^T
+                                 // \cdots B_{c'}^T B_c^T$ from the
+                                 // introduction that interpolates
+                                 // from the coarser of the two cells
+                                 // to the finer one. If the two cells
+                                 // match, then this is the identity
+                                 // matrix -- exactly what we pass to
+                                 // this function initially.
+                                 //
+                                 // The function has to consider two
+                                 // cases: that both of the two cells
+                                 // are not further refined, i.e. have
+                                 // no children, in which case we can
+                                 // finally assemble the right hand
+                                 // side contributions of this pair of
+                                 // cells; and that one of the two
+                                 // cells is further refined, in which
+                                 // case we have to keep recursing by
+                                 // looping over the children of the
+                                 // one cell that is not active. These
+                                 // two cases will be discussed below:
+template <int dim>
+void
+EnergyGroup<dim>::
+assemble_cross_group_rhs_recursive (const EnergyGroup<dim>                        &g_prime,
+                                   const typename DoFHandler<dim>::cell_iterator &cell_g,
+                                   const typename DoFHandler<dim>::cell_iterator &cell_g_prime,
+                                   const FullMatrix<double>                       prolongation_matrix)
+{
+                                   // The first case is that both
+                                   // cells are no further refined. In
+                                   // that case, we can assemble the
+                                   // relevant terms (see the
+                                   // introduction). This involves
+                                   // assembling the mass matrix on
+                                   // the finer of the two cells (in
+                                   // fact there are two mass matrices
+                                   // with different coefficients, one
+                                   // for the fission distribution
+                                   // cross section
+                                   // $\chi_g\nu\Sigma_{f,g'}$ and one
+                                   // for the scattering cross section
+                                   // $\Sigma_{s,g'\to g}$). This is
+                                   // straight forward, but note how
+                                   // we determine which of the two
+                                   // cells is ther finer one by
+                                   // looking at the refinement level
+                                   // of the two cells:
+  if (!cell_g->has_children() && !cell_g_prime->has_children())
+    {
+      const QGauss<dim>  quadrature_formula (fe.degree+1);
+      const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+    
+      FEValues<dim> fe_values (fe, quadrature_formula, 
+                               update_values  |  update_JxW_values);
+
+      if (cell_g->level() > cell_g_prime->level())
+       fe_values.reinit (cell_g);
+      else
+       fe_values.reinit (cell_g_prime);
+
+      const double fission_dist_XS
+       = material_data.get_fission_dist_XS (group, g_prime.group,
+                                             cell_g_prime->material_id());
+
+      const double scattering_XS
+       = material_data.get_scattering_XS (g_prime.group, group,
+                                           cell_g_prime->material_id());
+
+      FullMatrix<double>    local_mass_matrix_f (fe.dofs_per_cell,
+                                                 fe.dofs_per_cell);
+      FullMatrix<double>    local_mass_matrix_g (fe.dofs_per_cell,
+                                                 fe.dofs_per_cell);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+          for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+            {
+              local_mass_matrix_f(i,j) += (fission_dist_XS *
+                                           fe_values.shape_value(i,q_point) *
+                                           fe_values.shape_value(j,q_point) *
+                                           fe_values.JxW(q_point));
+              local_mass_matrix_g(i,j) += (scattering_XS *
+                                           fe_values.shape_value(i,q_point) *
+                                           fe_values.shape_value(j,q_point) *
+                                           fe_values.JxW(q_point));
+            }    
+
+                                       // Now we have all the
+                                       // interpolation (prolongation)
+                                       // matrices as well as local
+                                       // mass matrices, so we only
+                                       // have to form the product
+                                       // @f[
+                                       //  F_i|_{K_{cc'\cdots
+                                       //  c^{(k)}}} = [B_c B_{c'}
+                                       //  \cdots B_{c^{(k)}}
+                                       //  M_{K_{cc'\cdots
+                                       //  c^{(k)}}}]^{ij}
+                                       //  \phi_{g'}^j,
+                                       // @f]
+                                       // or
+                                       // @f[
+                                       //  F_i|_{K_{cc'\cdots
+                                       //  c^{(k)}}} = [(B_c B_{c'}
+                                       //  \cdots B_{c^{(k)}}
+                                       //  M_{K_{cc'\cdots
+                                       //  c^{(k)}}})^T]^{ij}
+                                       //  \phi_{g'}^j,
+                                       // @f]
+                                       // depending on which of the two
+                                       // cells is the finer. We do this
+                                       // using either the matrix-vector
+                                       // product provided by the <code>vmult</code>
+                                       // function, or the product with the
+                                       // transpose matrix using <code>Tvmult</code>.
+                                       // After doing so, we transfer the
+                                       // result into the global right hand
+                                       // side vector of energy group $g$.
+      Vector<double>       g_prime_new_values (fe.dofs_per_cell);
+      Vector<double>       g_prime_old_values (fe.dofs_per_cell);
+      cell_g_prime->get_dof_values (g_prime.solution_old, g_prime_old_values);
+      cell_g_prime->get_dof_values (g_prime.solution,     g_prime_new_values);
+      
+      Vector<double>       cell_rhs (fe.dofs_per_cell);
+      Vector<double>       tmp (fe.dofs_per_cell);
+
+      if (cell_g->level() > cell_g_prime->level())
+       {
+         prolongation_matrix.vmult (tmp, g_prime_old_values);
+         local_mass_matrix_f.vmult (cell_rhs, tmp);
+
+         prolongation_matrix.vmult (tmp, g_prime_new_values);
+         local_mass_matrix_g.vmult_add (cell_rhs, tmp);
+       }
+      else
+       {
+         local_mass_matrix_f.vmult (tmp, g_prime_old_values);
+         prolongation_matrix.Tvmult (cell_rhs, tmp);
+
+         local_mass_matrix_g.vmult (tmp, g_prime_new_values);
+         prolongation_matrix.Tvmult_add (cell_rhs, tmp);
+       }      
+
+      std::vector<unsigned int> local_dof_indices (fe.dofs_per_cell);
+      cell_g->get_dof_indices (local_dof_indices);
+      
+      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+       system_rhs(local_dof_indices[i]) += cell_rhs(i);
+    }
+
+                                   // The alternative is that one of
+                                   // the two cells is further
+                                   // refined. In that case, we have
+                                   // to loop over all the children,
+                                   // multiply the existing
+                                   // interpolation (prolongation)
+                                   // product of matrices from the
+                                   // left with the interpolation from
+                                   // the present cell to its child
+                                   // (using the matrix-matrix
+                                   // multiplication function
+                                   // <code>mmult</code>), and then hand the
+                                   // result off to this very same
+                                   // function again, but with the
+                                   // cell that has children replaced
+                                   // by one of its children:
+  else
+    for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell;++child)
+      {
+       FullMatrix<double>   new_matrix (fe.dofs_per_cell, fe.dofs_per_cell);
+       fe.get_prolongation_matrix(child).mmult (new_matrix,
+                                                 prolongation_matrix);
+
+       if (cell_g->has_children())
+         assemble_cross_group_rhs_recursive (g_prime,
+                                             cell_g->child(child), cell_g_prime,
+                                             new_matrix);
+       else
+         assemble_cross_group_rhs_recursive (g_prime,
+                                             cell_g, cell_g_prime->child(child),
+                                             new_matrix);
+      }
+}
+
+
+                                 // @sect5{<code>EnergyGroup::get_fission_source</code>}
+                                 //
+                                 // In the (inverse) power iteration,
+                                 // we use the integrated fission
+                                 // source to update the
+                                 // $k$-eigenvalue. Given its
+                                 // definition, the following function
+                                 // is essentially self-explanatory:
+template <int dim>
+double EnergyGroup<dim>::get_fission_source () const
+{
+  const QGauss<dim>  quadrature_formula (fe.degree + 1);
+  const unsigned int n_q_points    = quadrature_formula.n_quadrature_points;
+      
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values  |  update_JxW_values);
+      
+  std::vector<double>       solution_values (n_q_points);
+
+  double fission_source = 0;
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+
+      const double fission_XS
+       = material_data.get_fission_XS(group, cell->material_id());
+
+      fe_values.get_function_values (solution, solution_values);
+         
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       fission_source += (fission_XS *
+                          solution_values[q_point] *
+                          fe_values.JxW(q_point));
+    }
+      
+  return fission_source;
+}
+
+
+                                 // @sect5{<code>EnergyGroup::solve</code>}
+                                 //
+                                 // Next a function that solves the
+                                 // linear system assembled
+                                 // before. Things are pretty much
+                                 // standard, except that we delayed
+                                 // applying boundary values until we
+                                 // get here, since in all the
+                                 // previous functions we were still
+                                 // adding up contributions the right
+                                 // hand side vector.
+template <int dim>
+void
+EnergyGroup<dim>::solve () 
+{
+  hanging_node_constraints.condense (system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
+  
+  SolverControl           solver_control (system_matrix.m(),
+                                          1e-12*system_rhs.l2_norm());
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+  
+  cg.solve (system_matrix, solution, system_rhs, preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+}
+
+
+
+                                 // @sect5{<code>EnergyGroup::estimate_errors</code>}
+                                 //
+                                 // Mesh refinement is split into two
+                                 // functions. The first estimates the
+                                 // error for each cell, normalizes it
+                                 // by the magnitude of the solution,
+                                 // and returns it in the vector given
+                                 // as an argument. The calling
+                                 // function collects all error
+                                 // indicators from all energy groups,
+                                 // and computes thresholds for
+                                 // refining and coarsening cells.
+template <int dim>
+void EnergyGroup<dim>::estimate_errors (Vector<float> &error_indicators) const
+{
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+                                     QGauss<dim-1> (fe.degree + 1), 
+                                     typename FunctionMap<dim>::type(),
+                                     solution,
+                                     error_indicators);
+  error_indicators /= solution.linfty_norm();
+}
+
+
+
+                                 // @sect5{<code>EnergyGroup::refine_grid</code>}
+                                 //
+                                // The second part is to refine the
+                                // grid given the error indicators
+                                // compute in the previous function
+                                // and error thresholds above which
+                                // cells shall be refined or below
+                                // which cells shall be
+                                // coarsened. Note that we do not use
+                                // any of the functions in
+                                // <code>GridRefinement</code> here,
+                                // but rather set refinement flags
+                                // ourselves.
+                                //
+                                // After setting these flags, we use
+                                // the SolutionTransfer class to move
+                                // the solution vector from the old
+                                // to the new mesh. The procedure
+                                // used here is described in detail
+                                // in the documentation of that
+                                // class:
+template <int dim>
+void EnergyGroup<dim>::refine_grid (const Vector<float> &error_indicators,
+                                   const double         refine_threshold,
+                                   const double         coarsen_threshold)
+{
+  typename Triangulation<dim>::active_cell_iterator
+    cell = triangulation.begin_active(),
+    endc = triangulation.end();
+            
+  for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
+    if (error_indicators(cell_index) > refine_threshold)
+      cell->set_refine_flag ();
+    else if (error_indicators(cell_index) < coarsen_threshold)
+      cell->set_coarsen_flag ();
+      
+  SolutionTransfer<dim, double> soltrans(dof_handler);
+
+  triangulation.prepare_coarsening_and_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(solution);
+
+  triangulation.execute_coarsening_and_refinement ();
+  dof_handler.distribute_dofs (fe);
+
+  solution.reinit (dof_handler.n_dofs());
+  soltrans.interpolate(solution_old, solution);
+
+  solution_old.reinit (dof_handler.n_dofs());
+  solution_old = solution;
+}
+
+
+                                 // @sect5{<code>EnergyGroup::output_results</code>}
+                                 //
+                                 // The last function of this class
+                                 // outputs meshes and solutions after
+                                 // each mesh iteration. This has been
+                                 // shown many times before. The only
+                                 // thing worth pointing out is the
+                                 // use of the
+                                 // <code>Utilities::int_to_string</code>
+                                 // function to convert an integer
+                                 // into its string
+                                 // representation. The second
+                                 // argument of that function denotes
+                                 // how many digits we shall use -- if
+                                 // this value was larger than one,
+                                 // then the number would be padded by
+                                 // leading zeros.
+template <int dim>
+void
+EnergyGroup<dim>::output_results (const unsigned int cycle) const
+{
+  {
+    const std::string filename = std::string("grid-") +
+                                 Utilities::int_to_string(group,1) +
+                                 "." +
+                                 Utilities::int_to_string(cycle,1) +
+                                 ".eps";
+    std::ofstream output (filename.c_str());
+      
+    GridOut grid_out;
+    grid_out.write_eps (triangulation, output);
+  }
+
+  {
+    const std::string filename = std::string("solution-") +
+                                 Utilities::int_to_string(group,1) +
+                                 "." +
+                                 Utilities::int_to_string(cycle,1) +
+                                 ".gmv";
+       
+    DataOut<dim> data_out;
+
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, "solution");
+    data_out.build_patches ();
+  
+    std::ofstream output (filename.c_str());
+    data_out.write_gmv (output);
+  }
+}
+
+
+
+                                 // @sect3{The <code>NeutronDiffusionProblem</code> class template}
+
+                                 // This is the main class of the
+                                 // program, not because it implements
+                                 // all the functionality (in fact,
+                                 // most of it is implemented in the
+                                 // <code>EnergyGroup</code> class)
+                                 // but because it contains the
+                                 // driving algorithm that determines
+                                 // what to compute and when. It is
+                                 // mostly as shown in many of the
+                                 // other tutorial programs in that it
+                                 // has a public <code>run</code>
+                                 // function and private functions
+                                 // doing all the rest. In several
+                                 // places, we have to do something
+                                 // for all energy groups, in which
+                                 // case we will spawn threads for
+                                 // each group to let these things run
+                                 // in parallel if deal.II was
+                                 // configured for multithreading.
+                                 //
+                                 // The biggest difference to previous
+                                 // example programs is that we also
+                                 // declare a nested class that has
+                                 // member variables for all the
+                                 // run-time parameters that can be
+                                 // passed to the program in an input
+                                 // file. Right now, these are the
+                                 // number of energy groups, the
+                                 // number of refinement cycles, the
+                                 // polynomial degree of the finite
+                                 // element to be used, and the
+                                 // tolerance used to determine when
+                                 // convergence of the inverse power
+                                 // iteration has occurred. In
+                                 // addition, we have a constructor of
+                                 // this class that sets all these
+                                 // values to their default values, a
+                                 // function
+                                 // <code>declare_parameters</code>
+                                 // that described to the
+                                 // ParameterHandler class already
+                                 // used in @ref step_19 "step-19"
+                                 // what parameters are accepted in
+                                 // the input file, and a function
+                                 // <code>get_parameters</code> that
+                                 // can extract the values of these
+                                 // parameters from a ParameterHandler
+                                 // object.
+template <int dim>
+class NeutronDiffusionProblem 
+{
+  public:
+    class Parameters
+    {
+      public:
+       Parameters ();
+       
+       static void declare_parameters (ParameterHandler &prm);
+       void get_parameters (ParameterHandler &prm);
+
+       unsigned int n_groups;
+       unsigned int n_refinement_cycles;
+
+       unsigned int fe_degree;
+       
+       double convergence_tolerance;
+    };
+    
+
+    
+    NeutronDiffusionProblem (const Parameters &parameters);
+    ~NeutronDiffusionProblem ();
+
+    void run ();
+    
+  private:
+                                     // @sect5{Private member functions}
+    
+                                     // There are not that many member
+                                     // functions in this class since
+                                     // most of the functionality has
+                                     // been moved into the
+                                     // <code>EnergyGroup</code> class
+                                     // and is simply called from the
+                                     // <code>run()</code> member
+                                     // function of this class. The
+                                     // ones that remain have
+                                     // self-explanatory names:
+    void initialize_problem();
+
+    void refine_grid ();
+
+    double get_total_fission_source () const;
+
+
+                                     // @sect5{Private member variables}
+
+                                     // Next, we have a few member
+                                     // variables. In particular,
+                                     // these are (i) a reference to
+                                     // the parameter object (owned by
+                                     // the main function of this
+                                     // program, and passed to the
+                                     // constructor of this class),
+                                     // (ii) an object describing the
+                                     // material parameters for the
+                                     // number of energy groups
+                                     // requested in the input file,
+                                     // and (iii) the finite element
+                                     // to be used by all energy
+                                     // groups:
+    const Parameters  &parameters;
+    const MaterialData material_data;
+    FE_Q<dim>          fe;
+
+                                     // Furthermore, we have (iv) the
+                                     // value of the computed
+                                     // eigenvalue at the present
+                                     // iteration. This is, in fact,
+                                     // the only part of the solution
+                                     // that is shared between all
+                                     // energy groups -- all other
+                                     // parts of the solution, such as
+                                     // neutron fluxes are particular
+                                     // to one or the other energy
+                                     // group, and are therefore
+                                     // stored in objects that
+                                     // describe a single energy
+                                     // group:
+    double k_eff;
+
+                                     // Finally, (v), we have an array
+                                     // of pointers to the energy
+                                     // group objects. The length of
+                                     // this array is, of course,
+                                     // equal to the number of energy
+                                     // groups specified in the
+                                     // parameter file.
+    std::vector<EnergyGroup<dim>*> energy_groups;
+};
+
+
+                                 // @sect4{Implementation of the <code>NeutronDiffusionProblem::Parameters</code> class}
+
+                                // Before going on to the
+                                // implementation of the outer class,
+                                // we have to implement the functions
+                                // of the parameters structure. This
+                                // is pretty straightforward and, in
+                                // fact, looks pretty much the same
+                                // for all such parameters classes
+                                // using the ParameterHandler
+                                // capabilities. We will therefore
+                                // not comment further on this:
+template <int dim>
+NeutronDiffusionProblem<dim>::Parameters::Parameters ()
+               :
+               n_groups (2),
+               n_refinement_cycles (5),
+               fe_degree (2),
+               convergence_tolerance (1e-12)
+{}
+
+
+
+template <int dim>
+void
+NeutronDiffusionProblem<dim>::Parameters::
+declare_parameters (ParameterHandler &prm)
+{
+  prm.declare_entry ("Number of energy groups", "2",
+                    Patterns::Integer (),
+                    "The number of energy different groups considered");
+  prm.declare_entry ("Refinement cycles", "5",
+                    Patterns::Integer (),
+                    "Number of refinement cycles to be performed");
+  prm.declare_entry ("Finite element degree", "2",
+                    Patterns::Integer (),
+                    "Polynomial degree of the finite element to be used");
+  prm.declare_entry ("Power iteration tolerance", "1e-12",
+                    Patterns::Double (),
+                    "Inner power iterations are stopped when the change in k_eff falls "
+                    "below this tolerance");
+}
+
+
+
+template <int dim>
+void
+NeutronDiffusionProblem<dim>::Parameters::
+get_parameters (ParameterHandler &prm)
+{
+  n_groups              = prm.get_integer ("Number of energy groups");
+  n_refinement_cycles   = prm.get_integer ("Refinement cycles");
+  fe_degree             = prm.get_integer ("Finite element degree");
+  convergence_tolerance = prm.get_double ("Power iteration tolerance");
+}
+
+
+
+               
+                                 // @sect4{Implementation of the <code>NeutronDiffusionProblem</code> class}
+
+                                 // Now for the
+                                 // <code>NeutronDiffusionProblem</code>
+                                 // class. The constructor and
+                                 // destructor have nothing of much
+                                 // interest:
+template <int dim>
+NeutronDiffusionProblem<dim>::
+NeutronDiffusionProblem (const Parameters &parameters)
+               :
+               parameters (parameters),
+               material_data (parameters.n_groups),
+                fe (parameters.fe_degree)
+{}
+
+
+
+template <int dim>
+NeutronDiffusionProblem<dim>::~NeutronDiffusionProblem () 
+{
+  for (unsigned int group=0; group<energy_groups.size(); ++group)
+    delete energy_groups[group];
+
+  energy_groups.resize (0);
+}
+
+                                 // @sect5{<code>NeutronDiffusionProblem::initialize_problem</code>}
+                                 //
+                                 // The first function of interest is
+                                 // the one that sets up the geometry
+                                 // of the reactor core. This is
+                                 // described in more detail in the
+                                 // introduction.
+                                //
+                                // The first part of the function
+                                // defines geometry data, and then
+                                // creates a coarse mesh that has as
+                                // many cells as there are fuel rods
+                                // (or pin cells, for that matter) in
+                                // that part of the reactor core that
+                                // we simulate. As mentioned when
+                                // interpolating boundary values
+                                // above, the last parameter to the
+                                // <code>GridGenerator::subdivided_hyper_rectangle</code>
+                                // function specifies that sides of
+                                // the domain shall have unique
+                                // boundary indicators that will
+                                // later allow us to determine in a
+                                // simple way which of the boundaries
+                                // have Neumann and which have
+                                // Dirichlet conditions attached to
+                                // them.
+template <int dim>
+void NeutronDiffusionProblem<dim>::initialize_problem()
+{
+  const unsigned int rods_per_assembly_x = 17,
+                    rods_per_assembly_y = 17;
+  const double pin_pitch_x = 1.26,
+              pin_pitch_y = 1.26;
+  const double assembly_height = 200;
+
+  const unsigned int assemblies_x = 2,
+                    assemblies_y = 2,
+                    assemblies_z = 1;
+
+  const Point<dim> bottom_left = Point<dim>();
+  const Point<dim> upper_right = (dim == 2
+                                 ?
+                                 Point<dim> (assemblies_x*rods_per_assembly_x*pin_pitch_x,
+                                             assemblies_y*rods_per_assembly_y*pin_pitch_y)
+                                 :
+                                 Point<dim> (assemblies_x*rods_per_assembly_x*pin_pitch_x,
+                                             assemblies_y*rods_per_assembly_y*pin_pitch_y,
+                                             assemblies_z*assembly_height));     
+
+  std::vector<unsigned int> n_subdivisions;
+  n_subdivisions.push_back (assemblies_x*rods_per_assembly_x);
+  if (dim >= 2)
+    n_subdivisions.push_back (assemblies_y*rods_per_assembly_y);
+  if (dim >= 3)
+    n_subdivisions.push_back (assemblies_z);
+     
+  Triangulation<dim> coarse_grid;
+  GridGenerator::subdivided_hyper_rectangle (coarse_grid,
+                                             n_subdivisions,
+                                             bottom_left,
+                                             upper_right,
+                                             true);
+
+  
+                                   // The second part of the function
+                                   // deals with material numbers of
+                                   // pin cells of each type of
+                                   // assembly. Here, we define four
+                                   // different types of assembly, for
+                                   // which we describe the
+                                   // arrangement of fuel rods in the
+                                   // following tables.
+                                  //
+                                  // The assemblies described here
+                                  // are taken from the benchmark
+                                  // mentioned in the introduction
+                                  // and are (in this order):
+                                  // <ol>
+                                  //   <li>'UX' Assembly: UO2 fuel assembly
+                                  //       with 24 guide tubes and a central
+                                  //       Moveable Fission Chamber
+                                  //   <li>'UA' Assembly: UO2 fuel assembly
+                                  //       with 24 AIC and a central
+                                  //       Moveable Fission Chamber
+                                  //   <li>'PX' Assembly: MOX fuel assembly
+                                  //       with 24 guide tubes and a central
+                                  //       Moveable Fission Chamber
+                                  //   <li>'R' Assembly: a reflector.
+                                  // </ol>
+                                  //
+                                  // Note that the numbers listed
+                                  // here and taken from the
+                                  // benchmark description are, in
+                                  // good old Fortran fashion,
+                                  // one-based. We will later
+                                  // subtract one from each number
+                                  // when assigning materials to
+                                  // individual cells to convert
+                                  // things into the C-style
+                                  // zero-based indexing.
+  const unsigned int n_assemblies=4;
+  const unsigned int
+    assembly_materials[n_assemblies][rods_per_assembly_x][rods_per_assembly_y]
+    = {
+         {
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 5, 1, 1, 5, 1, 1, 7, 1, 1, 5, 1, 1, 5, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }
+         },      
+         {
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 8, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 8, 1, 1, 8, 1, 1, 7, 1, 1, 8, 1, 1, 8, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 8, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+               { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }
+         },
+         {
+               { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 },
+               { 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2 },
+               { 2, 3, 3, 3, 3, 5, 3, 3, 5, 3, 3, 5, 3, 3, 3, 3, 2 },
+               { 2, 3, 3, 5, 3, 4, 4, 4, 4, 4, 4, 4, 3, 5, 3, 3, 2 },
+               { 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2 },
+               { 2, 3, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 3, 2 },
+               { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+               { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+               { 2, 3, 5, 4, 4, 5, 4, 4, 7, 4, 4, 5, 4, 4, 5, 3, 2 },
+               { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+               { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+               { 2, 3, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 3, 2 },
+               { 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2 },
+               { 2, 3, 3, 5, 3, 4, 4, 4, 4, 4, 4, 4, 3, 5, 3, 3, 2 },
+               { 2, 3, 3, 3, 3, 5, 3, 3, 5, 3, 3, 5, 3, 3, 3, 3, 2 },
+               { 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2 },        
+               { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 }
+         },      
+         {
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+               { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 }
+         }
+    };
+
+                                  // After the description of the
+                                   // materials that make up an
+                                   // assembly, we have to specify the
+                                   // arrangement of assemblies within
+                                   // the core. We use a symmetric
+                                   // pattern that in fact only uses
+                                   // the 'UX' and 'PX' assemblies:
+  const unsigned int core[assemblies_x][assemblies_y][assemblies_z]
+    =  {{{0}, {2}}, {{2}, {0}}};
+  
+                                  // We are now in a position to
+                                  // actually set material IDs for
+                                  // each cell. To this end, we loop
+                                  // over all cells, look at the
+                                  // location of the cell's center,
+                                  // and determine which assembly and
+                                  // fuel rod this would be in. (We
+                                  // add a few checks to see that the
+                                  // locations we compute are within
+                                  // the bounds of the arrays in
+                                  // which we have to look up
+                                  // materials.) At the end of the
+                                  // loop, we set material
+                                  // identifiers accordingly:
+  for (typename Triangulation<dim>::active_cell_iterator
+         cell = coarse_grid.begin_active();
+       cell!=coarse_grid.end();
+       ++cell)
+    {
+      const Point<dim> cell_center = cell->center();
+        
+      const unsigned int tmp_x = int(cell_center[0]/pin_pitch_x);
+      const unsigned int ax = tmp_x/rods_per_assembly_x;
+      const unsigned int cx = tmp_x - ax * rods_per_assembly_x;
+
+      const unsigned tmp_y = int(cell_center[1]/pin_pitch_y);
+      const unsigned int ay = tmp_y/rods_per_assembly_y;
+      const unsigned int cy = tmp_y - ay * rods_per_assembly_y;
+
+      const unsigned int az = (dim == 2
+                              ?
+                              0
+                              :
+                              int (cell_center[dim-1]/assembly_height));
+
+      Assert (ax < assemblies_x, ExcInternalError());
+      Assert (ay < assemblies_y, ExcInternalError());
+      Assert (az < assemblies_z, ExcInternalError());
+
+      Assert (core[ax][ay][az] < n_assemblies, ExcInternalError());
+
+      Assert (cx < rods_per_assembly_x, ExcInternalError());
+      Assert (cy < rods_per_assembly_y, ExcInternalError());
+      
+      cell->set_material_id(assembly_materials[core[ax][ay][az]][cx][cy] - 1);
+    }
+
+                                  // With the coarse mesh so
+                                  // initialized, we create the
+                                  // appropriate number of energy
+                                  // group objects and let them
+                                  // initialize their individual
+                                  // meshes with the coarse mesh
+                                  // generated above:
+  energy_groups.resize (parameters.n_groups);
+  for (unsigned int group=0; group<parameters.n_groups; ++group)
+    energy_groups[group] = new EnergyGroup<dim> (group, material_data,
+                                                 coarse_grid, fe);
+}
+
+
+                                // @sect5{<code>NeutronDiffusionProblem::get_total_fission_source</code>}
+                                 //
+                                 // In the eigenvalue computation, we
+                                 // need to calculate total fission
+                                 // neutron source after each power
+                                 // iteration. The total power then is
+                                 // used to renew k-effective.
+                                //
+                                // Since the total fission source is
+                                // a sum over all the energy groups,
+                                // and since each of these sums can
+                                // be computed independently, we
+                                // actually do this in parallel. One
+                                // of the problems is that the
+                                // function in the
+                                // <code>EnergyGroup</code> class
+                                // that computes the fission source
+                                // returns a value. If we now simply
+                                // spin off a new thread, we have to
+                                // later capture the return value of
+                                // the function run on that
+                                // thread. The way this can be done
+                                // is to use the return value of the
+                                // Threads::spawn function, which is
+                                // of type Threads::Thread@<double@>
+                                // if the function spawned returns a
+                                // double. We can the later ask this
+                                // object for the returned value
+                                // (when doing so, the
+                                // Threads::Thread@<double@>::return_value
+                                // function first waits for the
+                                // thread to finish).
+                                //
+                                // The way this function then works
+                                // is to first spawn one thread for
+                                // each energy group we work with,
+                                // then one-by-one collecting the
+                                // returned values of each thread and
+                                // return the sum.
+template <int dim>
+double NeutronDiffusionProblem<dim>::get_total_fission_source () const
+{
+  std::vector<Threads::Thread<double> > threads;
+  for (unsigned int group=0; group<parameters.n_groups; ++group)
+    threads.push_back (Threads::spawn (*energy_groups[group],
+                                      &EnergyGroup<dim>::get_fission_source) ());
+  
+  double fission_source = 0;
+  for (unsigned int group=0; group<parameters.n_groups; ++group)
+    fission_source += threads[group].return_value ();
+
+  return fission_source;
+}
+
+
+
+
+                                // @sect5{<code>NeutronDiffusionProblem::refine_grid</code>}
+                                 //
+                                 // The next function lets the
+                                 // individual energy group objects
+                                 // refine their meshes. Much of this,
+                                 // again, is a task that can be done
+                                 // independently in parallel: first,
+                                 // let all the energy group objects
+                                 // calculate their error indicators
+                                 // in parallel, then compute the
+                                 // maximum error indicator over all
+                                 // energy groups and determine
+                                 // thresholds for refinement and
+                                 // coarsening of cells, and then ask
+                                 // all the energy groups to refine
+                                 // their meshes accordingly, again in
+                                 // parallel.
+template <int dim>
+void NeutronDiffusionProblem<dim>::refine_grid ()
+{
+  std::vector<unsigned int> n_cells (parameters.n_groups);
+  for (unsigned int group=0; group<parameters.n_groups; ++group)
+    n_cells[group] = energy_groups[group]->n_active_cells();
+
+  BlockVector<float>  group_error_indicators(n_cells);
+
+  {  
+    Threads::ThreadGroup<> threads;
+    for (unsigned int group=0; group<parameters.n_groups; ++group)
+      threads += Threads::spawn (*energy_groups[group], &EnergyGroup<dim>::estimate_errors)
+                (group_error_indicators.block(group));
+    threads.join_all ();
+  }
+  
+  const float max_error         = group_error_indicators.linfty_norm();
+  const float refine_threshold  = 0.3*max_error;
+  const float coarsen_threshold = 0.01*max_error;
+
+  {  
+    Threads::ThreadGroup<> threads;
+    for (unsigned int group=0; group<parameters.n_groups; ++group)
+      threads += Threads::spawn (*energy_groups[group], &EnergyGroup<dim>::refine_grid)
+                (group_error_indicators.block(group),
+                 refine_threshold,
+                 coarsen_threshold);
+    threads.join_all ();
+  }
+}
+
+
+                                // @sect5{<code>NeutronDiffusionProblem::run</code>}
+                                 //
+                                 // Finally, this is the function
+                                 // where the meat is: iterate on a
+                                 // sequence of meshes, and on each of
+                                 // them do a power iteration to
+                                 // compute the eigenvalue.
+                                //
+                                // Given the description of the
+                                // algorithm in the introduction,
+                                // there is actually not much to
+                                // comment on:
+template <int dim>
+void NeutronDiffusionProblem<dim>::run () 
+{
+  std::cout << std::setprecision (12) << std::fixed;
+
+  double k_eff_old = k_eff;
+
+  Timer timer;
+  timer.start ();
+  
+  for (unsigned int cycle=0; cycle<parameters.n_refinement_cycles; ++cycle)
+    {
+      std::cout << "Cycle " << cycle << ':' << std::endl;
+      
+      if (cycle == 0)
+       initialize_problem();
+      else
+       {
+         refine_grid ();
+         for (unsigned int group=0; group<parameters.n_groups; ++group)
+           energy_groups[group]->solution *= k_eff;
+       } 
+
+      for (unsigned int group=0; group<parameters.n_groups; ++group)
+       energy_groups[group]->setup_linear_system ();
+
+      std::cout << "   Numbers of active cells:       ";
+      for (unsigned int group=0; group<parameters.n_groups; ++group)
+       std::cout << energy_groups[group]->n_active_cells()
+                 << ' ';
+      std::cout << std::endl;
+      std::cout << "   Numbers of degrees of freedom: ";
+      for (unsigned int group=0; group<parameters.n_groups; ++group)
+       std::cout << energy_groups[group]->n_dofs()
+                 << ' ';
+      std::cout << std::endl << std::endl;
+      
+      
+      Threads::ThreadGroup<> threads;
+      for (unsigned int group=0; group<parameters.n_groups; ++group)
+       threads += Threads::spawn
+                  (*energy_groups[group], &EnergyGroup<dim>::assemble_system_matrix)
+                  ();
+      threads.join_all ();
+
+      double max_old = 0;
+
+                                       // indicate this is a eigenvalue problem
+      unsigned int isour = 0;
+                                       // store relative error between two
+                                       // successive power iterations
+      double error;
+
+      unsigned int iteration = 1;
+      do
+       {
+         for (unsigned int group=0; group<parameters.n_groups; ++group)
+           {
+             energy_groups[group]->assemble_ingroup_rhs (ZeroFunction<dim>());
+             
+             for (unsigned int bgroup=0; bgroup<parameters.n_groups; ++bgroup)
+               energy_groups[group]->assemble_cross_group_rhs (*energy_groups[bgroup]);
+             
+             energy_groups[group]->solve ();
+           }
+
+         if (isour==1)
+           {
+             double max_current = 0;
+              for (unsigned int group=0; group<parameters.n_groups; ++group)
+               {             
+                 max_current = std::max (max_current,
+                                         energy_groups[group]->solution.linfty_norm());
+                 
+                 energy_groups[group]->solution_old = energy_groups[group]->solution;
+               }
+
+             error = fabs(max_current-max_old)/max_current;
+             max_old = max_current;
+           }
+         else
+           { 
+             k_eff = get_total_fission_source();
+             error = fabs(k_eff-k_eff_old)/fabs(k_eff);
+             std::cout << "   Iteration " << iteration
+                       << ": k_eff=" << k_eff
+                       << std::endl;
+             k_eff_old=k_eff;
+
+              for (unsigned int group=0; group<parameters.n_groups; ++group)
+                {
+                  energy_groups[group]->solution_old = energy_groups[group]->solution;
+                  energy_groups[group]->solution_old /= k_eff;
+                }
+           }
+
+         ++iteration;
+       }
+      while((error > parameters.convergence_tolerance)
+           &&
+           (iteration < 500));
+             
+      for (unsigned int group=0; group<parameters.n_groups; ++group)
+        energy_groups[group]->output_results (cycle);
+
+      std::cout << std::endl;
+      std::cout << "   Cycle=" << cycle
+               << ", n_dofs=" << energy_groups[0]->n_dofs() + energy_groups[1]->n_dofs()
+               << ",  k_eff=" << k_eff
+               << ", time=" << timer()
+               << std::endl;
+
+      
+      std::cout << std::endl << std::endl;
+    }
+}
+
+
+
+                                // @sect3{The <code>main()</code> function}
+                                 //
+                                // The last thing in the program in
+                                // the <code>main()</code>
+                                // function. The structure is as in
+                                // most other tutorial programs, with
+                                // the only exception that we here
+                                // handle a parameter file.  To this
+                                // end, we first look at the command
+                                // line arguments passed to this
+                                // function: if no input file is
+                                // specified on the command line,
+                                // then use "project.prm", otherwise
+                                // take the filename given as the
+                                // first argument on the command
+                                // line.
+                                //
+                                // With this, we create a
+                                // ParameterHandler object, let the
+                                // <code>NeutronDiffusionProblem::Parameters</code>
+                                // class declare all the parameters
+                                // it wants to see in the input file
+                                // (or, take the default values, if
+                                // nothing is listed in the parameter
+                                // file), then read the input file,
+                                // ask the parameters object to
+                                // extract the values, and finally
+                                // hand everything off to an object
+                                // of type
+                                // <code>NeutronDiffusionProblem</code>
+                                // for computation of the eigenvalue:
+int main (int argc, char ** argv) 
+{
+  const unsigned int dim = 2;
+  
+  try
+    {
+      deallog.depth_console (0);
+
+      std::string filename;
+      if (argc < 2)
+       filename = "project.prm";
+      else
+       filename = argv[1];
+
+      
+      ParameterHandler parameter_handler;
+      
+      NeutronDiffusionProblem<dim>::Parameters parameters;
+      parameters.declare_parameters (parameter_handler);
+
+      parameter_handler.read_input (filename);
+      
+      parameters.get_parameters (parameter_handler);
+
+      
+      NeutronDiffusionProblem<dim> neutron_diffusion_problem (parameters);
+      neutron_diffusion_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;
+}
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.