]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add step-23, 24, and 25 in its present, partly undocumented, form
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Jul 2006 16:13:04 +0000 (16:13 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Jul 2006 16:13:04 +0000 (16:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@13431 0785d39b-7218-0410-832d-ea1e28bc413d

37 files changed:
deal.II/examples/step-23/Makefile [new file with mode: 0644]
deal.II/examples/step-23/Makefile.dep [new file with mode: 0644]
deal.II/examples/step-23/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-23/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-23/step-23 [new file with mode: 0755]
deal.II/examples/step-23/step-23.cc [new file with mode: 0644]
deal.II/examples/step-24/Makefile [new file with mode: 0644]
deal.II/examples/step-24/Makefile.dep [new file with mode: 0644]
deal.II/examples/step-24/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-24/doc/multi.png [new file with mode: 0644]
deal.II/examples/step-24/doc/multi_movie.gif [new file with mode: 0644]
deal.II/examples/step-24/doc/multi_s.png [new file with mode: 0644]
deal.II/examples/step-24/doc/multi_s2.png [new file with mode: 0644]
deal.II/examples/step-24/doc/multi_sf.png [new file with mode: 0644]
deal.II/examples/step-24/doc/one.png [new file with mode: 0644]
deal.II/examples/step-24/doc/one_movie.gif [new file with mode: 0644]
deal.II/examples/step-24/doc/one_s.png [new file with mode: 0644]
deal.II/examples/step-24/doc/one_s2.png [new file with mode: 0644]
deal.II/examples/step-24/doc/one_sf.png [new file with mode: 0644]
deal.II/examples/step-24/doc/project-1.tex [new file with mode: 0644]
deal.II/examples/step-24/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-24/step-24 [new file with mode: 0755]
deal.II/examples/step-24/step-24.cc [new file with mode: 0644]
deal.II/examples/step-25/Makefile [new file with mode: 0644]
deal.II/examples/step-25/Makefile.dep [new file with mode: 0644]
deal.II/examples/step-25/doc/animate.sh [new file with mode: 0755]
deal.II/examples/step-25/doc/christov.1d-breather.gif [new file with mode: 0644]
deal.II/examples/step-25/doc/christov.1d-breather_stp.png [new file with mode: 0644]
deal.II/examples/step-25/doc/christov.2d-angled_kink.gif [new file with mode: 0644]
deal.II/examples/step-25/doc/christov.2d-kink.gif [new file with mode: 0644]
deal.II/examples/step-25/doc/christov.2d-kink.png [new file with mode: 0644]
deal.II/examples/step-25/doc/christov.2d-pseudobreather.gif [new file with mode: 0644]
deal.II/examples/step-25/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-25/doc/plot.plt [new file with mode: 0644]
deal.II/examples/step-25/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-25/step-25 [new file with mode: 0755]
deal.II/examples/step-25/step-25.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-23/Makefile b/deal.II/examples/step-23/Makefile
new file mode 100644 (file)
index 0000000..f291d79
--- /dev/null
@@ -0,0 +1,154 @@
+# $Id: Makefile 11909 2005-12-21 13:30:36Z guido $
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+################################################################
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<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-23/Makefile.dep b/deal.II/examples/step-23/Makefile.dep
new file mode 100644 (file)
index 0000000..f6e3599
--- /dev/null
@@ -0,0 +1,242 @@
+./step-23.o: \
+               step-23.cc\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/config.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/data_out_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function_time.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/geometry_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/logstream.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/multithread_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/numbers.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/point.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/polynomial.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/smartpointer.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/subscriptor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_product_polynomials.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/thread_management.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/vector_slice.h\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/assert.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/checked_delete.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/posix_features.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_compiler_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_platform_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_stdlib_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/suffix.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/current_function.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_gcc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_win32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/bad_weak_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/interlocked.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lightweight_mutex.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_nop.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_win32_cs.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/quick_allocator.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_ptr_nmt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_cw_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ia64.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_x86.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_nt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_pt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_w32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_impl.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/workaround.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/shared_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/throw_exception.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/tuple/tuple.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/alignment_of.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/type_with_alignment.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_constraints.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/function_map.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_collection.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_poly.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_q.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_update_flags.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_values.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_generator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_faces.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_hex.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_line.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_quad.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/matrices.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/vectors.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/full_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/identity_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/precondition.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_cg.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_control.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparse_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparsity_pattern.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector_memory.h
+./step-23.g.o: \
+               step-23.cc\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/config.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/data_out_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function_time.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/geometry_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/logstream.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/multithread_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/numbers.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/point.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/polynomial.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/smartpointer.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/subscriptor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_product_polynomials.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/thread_management.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/vector_slice.h\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/assert.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/checked_delete.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/posix_features.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_compiler_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_platform_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_stdlib_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/suffix.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/current_function.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_gcc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_win32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/bad_weak_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/interlocked.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lightweight_mutex.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_nop.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_win32_cs.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/quick_allocator.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_ptr_nmt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_cw_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ia64.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_x86.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_nt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_pt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_w32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_impl.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/workaround.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/shared_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/throw_exception.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/tuple/tuple.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/alignment_of.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/type_with_alignment.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_constraints.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/function_map.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_collection.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_poly.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_q.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_update_flags.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_values.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_generator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_faces.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_hex.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_line.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_quad.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/matrices.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/vectors.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/full_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/identity_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/precondition.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_cg.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_control.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparse_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparsity_pattern.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector_memory.h
diff --git a/deal.II/examples/step-23/doc/intro.dox b/deal.II/examples/step-23/doc/intro.dox
new file mode 100644 (file)
index 0000000..c59cfdf
--- /dev/null
@@ -0,0 +1,2 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
diff --git a/deal.II/examples/step-23/doc/results.dox b/deal.II/examples/step-23/doc/results.dox
new file mode 100644 (file)
index 0000000..e67fccc
--- /dev/null
@@ -0,0 +1,2 @@
+<a name="Results"></a>
+<h1>Results</h1>
diff --git a/deal.II/examples/step-23/step-23 b/deal.II/examples/step-23/step-23
new file mode 100755 (executable)
index 0000000..21366a8
Binary files /dev/null and b/deal.II/examples/step-23/step-23 differ
diff --git a/deal.II/examples/step-23/step-23.cc b/deal.II/examples/step-23/step-23.cc
new file mode 100644 (file)
index 0000000..65ffc0c
--- /dev/null
@@ -0,0 +1,357 @@
+/* $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */
+
+/*    $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $       */
+/*    Version: $Name:  $                                          */
+/*                                                                */
+/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors */
+/*                                                                */
+/*    This file is subject to QPL and may not be  distributed     */
+/*    without copyright and license information. Please refer     */
+/*    to the file deal.II/doc/license.html for the  text  and     */
+/*    further information on this license.                        */
+
+
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <dofs/dof_constraints.h>
+
+#include <numerics/matrices.h>
+#include <numerics/vectors.h>
+
+#include <numerics/data_out.h>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#include <base/logstream.h>
+
+
+
+template <int dim>
+class WaveEquationProblem 
+{
+  public:
+    WaveEquationProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve_u ();
+    void solve_v ();
+    void output_results (const unsigned int timestep_number) const;
+
+    Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    SparseMatrix<double> mass_matrix;
+    SparseMatrix<double> laplace_matrix;
+
+    double time_step;
+    double theta;
+
+    Vector<double>       solution_u, solution_v;
+    Vector<double>       old_solution_u, old_solution_v;
+    Vector<double>       system_rhs;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim> 
+{
+  public:
+    RightHandSide () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+template <int dim>
+class InitialValuesU : public Function<dim> 
+{
+  public:
+    InitialValuesU () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim> 
+{
+  public:
+    BoundaryValues () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &/*p*/,
+                                 const unsigned int /*component*/) const 
+{
+//   if (get_time() <= 0.25)
+//     if ((p[0] <=0) && (p[1] <= 0))
+//       return 1;
+
+  return 0;
+}
+
+
+template <int dim>
+double InitialValuesU<dim>::value (const Point<dim> &p,
+                                 const unsigned int /*component*/) const 
+{
+  //  return std::sqrt(p.square()) * std::exp (-p.square()) / 3;
+  if ((p[0] <=0) && (p[1] <= 0))
+    return 1;
+
+  return 0;}
+
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &/*p*/,
+                                  const unsigned int /*component*/) const 
+{
+  return 0;
+}
+
+
+
+
+
+
+
+template <int dim>
+WaveEquationProblem<dim>::WaveEquationProblem () :
+                fe (1),
+               dof_handler (triangulation),
+               time_step (1./64),
+               theta (0.5)
+{}
+
+
+
+template <int dim>
+void WaveEquationProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (7);
+  
+  std::cout << "   Number of active cells: "
+           << triangulation.n_active_cells()
+           << std::endl
+           << "   Total number of cells: "
+           << triangulation.n_cells()
+           << std::endl;
+
+  dof_handler.distribute_dofs (fe);
+
+  std::cout << "   Number of degrees of freedom: "
+           << dof_handler.n_dofs()
+           << std::endl;
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+  mass_matrix.reinit (sparsity_pattern);
+  laplace_matrix.reinit (sparsity_pattern);
+
+  solution_u.reinit (dof_handler.n_dofs());
+  solution_v.reinit (dof_handler.n_dofs());
+  old_solution_u.reinit (dof_handler.n_dofs());
+  old_solution_v.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void WaveEquationProblem<dim>::assemble_system () 
+{  
+  MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3),
+                                    mass_matrix);
+  MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3),
+                                       laplace_matrix);
+  
+  system_matrix.copy_from (mass_matrix);
+  system_matrix.add (theta * theta * time_step * time_step, laplace_matrix);
+}
+
+
+
+template <int dim>
+void WaveEquationProblem<dim>::solve_u () 
+{
+  SolverControl           solver_control (1000, 1e-8*system_rhs.l2_norm());
+  SolverCG<>              cg (solver_control);
+  cg.solve (system_matrix, solution_u, system_rhs,
+           PreconditionIdentity());
+
+  std::cout << "   " << solver_control.last_step()
+           << " CG iterations needed to obtain convergence."
+           << std::endl;
+}
+
+
+template <int dim>
+void WaveEquationProblem<dim>::solve_v () 
+{
+  SolverControl           solver_control (1000, 1e-8*system_rhs.l2_norm());
+  SolverCG<>              cg (solver_control);
+  cg.solve (mass_matrix, solution_v, system_rhs,
+           PreconditionIdentity());
+
+  std::cout << "   " << solver_control.last_step()
+           << " CG iterations needed to obtain convergence."
+           << std::endl;
+}
+
+
+
+template <int dim>
+void WaveEquationProblem<dim>::output_results (const unsigned int timestep_number) const
+{
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution_u, "U");
+  data_out.add_data_vector (solution_v, "V");
+
+  data_out.build_patches ();
+
+  std::ostringstream filename;
+  filename << "solution-"
+          << timestep_number;
+  std::ofstream output (filename.str().c_str());
+  data_out.write_gnuplot (output);
+}
+
+
+
+
+template <int dim>
+void WaveEquationProblem<dim>::run () 
+{
+  std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
+  
+  make_grid_and_dofs();
+  assemble_system ();
+
+  ConstraintMatrix constraints;
+  constraints.close();
+  VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
+                       InitialValuesU<dim>(),
+                       old_solution_u);
+  VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
+                       ZeroFunction<dim>(),
+                       old_solution_v);
+
+  unsigned int timestep_number = 1;
+  for (double time = time_step; time<=5; time+=time_step, ++timestep_number)
+    {
+      Vector<double> tmp (solution_u.size());
+
+      mass_matrix.vmult (system_rhs, old_solution_u);
+
+      mass_matrix.vmult (tmp, old_solution_v);
+      system_rhs.add (time_step, tmp);
+
+      laplace_matrix.vmult (tmp, old_solution_u);
+      system_rhs.add (-theta * (1-theta) * time_step * time_step, tmp);
+
+      RightHandSide<dim> rhs_function;
+      rhs_function.set_time (time);
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                          rhs_function, tmp);
+      system_rhs.add (theta * theta * time_step * time_step, tmp);
+
+      rhs_function.set_time (time-time_step);
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                          rhs_function, tmp);
+      system_rhs.add (theta * (1-theta) * time_step * time_step, tmp);
+
+
+      std::map<unsigned int,double> boundary_values;
+      VectorTools::interpolate_boundary_values (dof_handler,
+                                               0,
+                                               BoundaryValues<dim>(),
+                                               boundary_values);
+      MatrixTools::apply_boundary_values (boundary_values,
+                                         system_matrix,
+                                         solution_u,
+                                         system_rhs);
+      solve_u ();
+
+
+      laplace_matrix.vmult (system_rhs, solution_u);
+      system_rhs *= -theta * time_step;
+
+      mass_matrix.vmult (tmp, old_solution_v);
+      system_rhs += tmp;
+
+      laplace_matrix.vmult (tmp, old_solution_u);
+      system_rhs.add (-time_step * (1-theta), tmp);
+
+      rhs_function.set_time (time);
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                          rhs_function, tmp);
+      system_rhs.add (theta * time_step, tmp);
+
+      rhs_function.set_time (time-time_step);
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                          rhs_function, tmp);
+      system_rhs.add ((1-theta) * time_step, tmp);
+
+      solve_v ();
+
+      output_results (timestep_number);
+
+      old_solution_u = solution_u;
+      old_solution_v = solution_v;
+    }
+}
+
+
+
+int main () 
+{
+  deallog.depth_console (0);
+  {
+    WaveEquationProblem<2> wave_equation_problem_2d;
+    wave_equation_problem_2d.run ();
+  }
+  
+  return 0;
+}
diff --git a/deal.II/examples/step-24/Makefile b/deal.II/examples/step-24/Makefile
new file mode 100644 (file)
index 0000000..f291d79
--- /dev/null
@@ -0,0 +1,154 @@
+# $Id: Makefile 11909 2005-12-21 13:30:36Z guido $
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+################################################################
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<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-24/Makefile.dep b/deal.II/examples/step-24/Makefile.dep
new file mode 100644 (file)
index 0000000..94cc50a
--- /dev/null
@@ -0,0 +1,252 @@
+./step-24.o: \
+               step-24.cc\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/config.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/data_out_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function_time.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/geometry_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/logstream.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/multithread_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/numbers.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/point.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/polynomial.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/qprojector.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/smartpointer.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/subscriptor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_product_polynomials.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/thread_management.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/vector_slice.h\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/assert.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/checked_delete.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/posix_features.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_compiler_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_platform_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_stdlib_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/suffix.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/current_function.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_gcc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_win32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/bad_weak_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/interlocked.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lightweight_mutex.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_nop.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_win32_cs.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/quick_allocator.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_ptr_nmt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_cw_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ia64.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_x86.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_nt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_pt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_w32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_impl.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/workaround.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/shared_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/throw_exception.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/tuple/tuple.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/alignment_of.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/type_with_alignment.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_constraints.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/function_map.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_collection.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_poly.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_q.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_update_flags.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_values.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping_q1.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_generator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_boundary.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_boundary_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_faces.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_hex.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_line.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_quad.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/matrices.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/vectors.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/full_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/identity_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/precondition.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_cg.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_control.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparse_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparsity_pattern.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector_memory.h
+./step-24.g.o: \
+               step-24.cc\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/config.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/data_out_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function_time.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/geometry_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/logstream.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/multithread_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/numbers.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/point.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/polynomial.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/qprojector.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/smartpointer.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/subscriptor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_product_polynomials.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/thread_management.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/vector_slice.h\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/assert.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/checked_delete.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/posix_features.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_compiler_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_platform_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_stdlib_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/suffix.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/current_function.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_gcc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_win32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/bad_weak_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/interlocked.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lightweight_mutex.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_nop.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_win32_cs.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/quick_allocator.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_ptr_nmt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_cw_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ia64.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_x86.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_nt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_pt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_w32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_impl.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/workaround.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/shared_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/throw_exception.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/tuple/tuple.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/alignment_of.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/type_with_alignment.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_constraints.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/function_map.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_collection.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_poly.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_q.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_update_flags.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_values.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping_q1.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_generator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_boundary.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_boundary_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_faces.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_hex.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_line.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_quad.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/matrices.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/vectors.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/full_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/identity_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/precondition.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_cg.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_control.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparse_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparsity_pattern.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector_memory.h
diff --git a/deal.II/examples/step-24/doc/intro.dox b/deal.II/examples/step-24/doc/intro.dox
new file mode 100644 (file)
index 0000000..939112b
--- /dev/null
@@ -0,0 +1,162 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+This project is to simulate the thermoacoustic tomography imaging. In thermoacoustic
+tomography,pulsed electromagnetic energy is delivered into biological issues. 
+Tissues absorbe the energy and then generate thermoacoustic waves through 
+thermoelastic expansion. The forward problem is a wave propagation problem.
+
+
+<h3>Problem</h3>
+
+Thermal equation without considering thermal diffusion is
+
+@f[
+\rho C_p \frac{\partial}{\partial t}T(t,\mathbf r) = H(t,\mathbf r)
+@f]
+
+Where $\rho (\mathbf r) $ is the density; $C_p (\mathbf r) $ is the specific heat; 
+$ T(t,\mathbf r)$ is the temperature rise due to the delivered microwave
+energy; and $H(t,\mathbf r)$ is the heating function defined as the thermal energy 
+per time and volume transformed from deposited microwave energy.
+
+Assume tissues have heterogeneous dielectric properties but homogeneous acoustic 
+properties. The basic acoustic generation equation in an acoustically homogeneous
+medium is the linear inviscid force equation
+
+@f[
+\rho \frac{\partial^2}{\partial t^2}u(t,\mathbf r) = 
+-\nabla p(t,\mathbf r)
+@f]
+
+and the expansion equation:
+
+@f[
+\nabla \cdot u(t,\mathbf r) = -\frac{p(t,\mathbf r)}{\rho c_0^2}+\beta T(t,\mathbf r) 
+@f]
+
+The original problem can be described as:
+
+@f[
+\Delta p-\frac{1}{c_0^2} \frac{\partial^2 p   }{\partial^2 t} = \lambda \delta(t)a(\mathbf r)
+@f]
+
+where $\lambda = - \frac{\beta}{C_p}$.
+
+The forward propogation problem can be changed to solve a wave equation with 
+initial conditions as follows:
+
+@f{eqnarray*} 
+\Delta \bar{p}- \frac{1}{c_0^2} \frac{\partial^2 \bar{p}}{\partial^2 t} & = & f(t,\mathbf r) \\
+
+\bar{p}(0,\mathbf r)=\lambda a(\mathbf r) & = & b(\mathbf r) 
+@f}
+
+
+
+<h3>Weak form and Discretization</h3>
+
+One first introduces a second variable, which is defined as the derivative of
+the pressure potential.
+
+@f[ 
+v = \frac{\partial\bar{p}}{\partial t} 
+@f]
+
+With the second variables, one then transform the forward problem into
+two seperate equations:
+
+@f{eqnarray*}
+\bar{p}_{t} - v & = & 0 \\
+\Delta\bar{p} - \frac{1}{c_0^2}\,v_{t} & = & f 
+@f}
+
+with initial conditions:
+
+@f{eqnarray*}
+\bar{p}(0,\mathbf r) & = & b(r) \\
+v(0,\mathbf r)=\bar{p}_t(0,\mathbf r) & = & 1 
+@f}
+
+In real application, the thermoacoustic source is very small as compared to the medium.  
+The propagation path of the thermoacoustic waves can be approximated as from the source 
+to the infinity. And the detector is in limited distance from the source. One only needs to
+evaluate the values when the thermoacoustic waves pass through the detectors. For this specific
+detection geometry, One then chooses the absorbing boundary condition for the simulation.
+
+@f[
+\frac{\partial\bar{p}}{\partial\mathbf n} = 
+-\frac{1}{c_0} \frac{\partial\bar{p}}{\partial t}
+@f]
+
+$\frac{\partial\bar{p}}{\partial\mathbf n}$ is normal derivative at the boundary. This is the 
+time-varying FEM model. To implement FEM for time dependent problem , one discretizes according    
+to $t$ and obtains: 
+
+@f{eqnarray*}
+(\frac{\bar{p}^n-\bar{p}^{n-1}}{\delta t},\phi)_\Omega-
+(\theta v^{n}+(1-\theta)v^{n-1},\phi)_\Omega & = & 0   \\
+-(\Delta((\theta\bar{p}^n+(1-\theta)\bar{p}^{n-1})),\nabla)_\Omega-
+\frac{1}{c_0}(\frac{\bar{p}^n-\bar{p}^{n-1}}{\delta t},\phi)_\partial\Omega - 
+\frac{1}{c_0^2}(\frac{v^n-v^{n-1}}{\delta t},\phi)_\Omega & = 
+& \theta f^{n}+(1-\theta)f^{n-1}
+@f}
+
+The weak formulation of the problem is obtained by multiplying the above two equations
+with test functions and integrating some terms by parts:
+
+@f{eqnarray*}
+M\bar{p}^{n}-(\delta t \theta)M v^{n-1} & = & M\bar{p}^{n-1}+\delta t (1-\theta)Mv^{n-1}\\
+
+(-c_0^2\delta t \theta A-c_0 B)\bar{p}^n-Mv^{n} & = & 
+(c_0^2\delta t(1-\theta)A-c_0B)\bar{p}^{n-1}-Mv^{n-1}+c_0^2\delta t(\theta F^{n}+(1-\theta)F^{n-1})
+@f}
+
+Here, the absoring boundary conditions are incorporated into the weak form by using 
+
+@f[ 
+\int_\Omega\varphi(\nabla\cdot(\nabla p))dx =
+-\int_\Omega\nabla \varphi \cdot \nabla p dx + 
+\int_{\partial\Omega}\varphi \frac{\partial p}{\partial t}ds
+@f]
+
+Where $\varphi$ is the test function.
+
+Pressure and its derivative are the two variables one is interested in the above equations,
+One can write the above two equations as a matrix form with the pressure and its derivative as
+an unknown vecotor:
+@f[
+\left(\begin{array}{cc}
+ M         &       -\delta t\theta M \\
+c_0^2\,\delta t\,\theta\,A+c_0\,B  &  M   \\
+               \end{array} \right)\\
+\left(\begin{array}{c}
+ \bar{p}^{n}    \\
+ \bar{v}^{n}
+              \end{array}\right)=\\
+\left(\begin{array}{l}
+ G_1  \\
+ G_2 -(\theta F^{n}+(1-\theta)F ^{n-1})c_{0}^{2}\delta t \\        
+                \end{array}\right)
+@f]
+
+where
+@f[
+\left(\begin{array}{c}
+G_1 \\
+G_2 \\
+   \end{array} \right)=\\
+\left(\begin{array}{l}
+ M\bar{p}^{n-1}+\delta t(1-\theta)Mv^{n-1}\\               
+ (-c_{0}^{2}\delta t (1-\theta)A+c_0 B)\bar{p}^{n-1} +Mv^{n-1}
+                \end{array}\right)
+@f]
+
+By some simply transformation, one obtains two iterative equations for
+the pressure potential and its derivative:
+@f{eqnarray*}
+(M+(\delta t\,\theta\,c_{0})^{2}A+c_0\delta t\theta B)\bar{p}^{n} & = & 
+G_{1}+(\delta t\, \theta)G_{2}-(c_0\delta t)^2\theta (\theta F^{n}+(1-\theta)F^{n-1}) \\
+Mv^n & = & -(c_0^2\,\delta t\, \theta\, A+c_0B)\bar{p}^{n}+ G_2 - 
+c_0^2\delta t(\theta F^{n}+(1-\theta)F^{n-1}) 
+@f}
diff --git a/deal.II/examples/step-24/doc/multi.png b/deal.II/examples/step-24/doc/multi.png
new file mode 100644 (file)
index 0000000..509201d
Binary files /dev/null and b/deal.II/examples/step-24/doc/multi.png differ
diff --git a/deal.II/examples/step-24/doc/multi_movie.gif b/deal.II/examples/step-24/doc/multi_movie.gif
new file mode 100644 (file)
index 0000000..c922c67
Binary files /dev/null and b/deal.II/examples/step-24/doc/multi_movie.gif differ
diff --git a/deal.II/examples/step-24/doc/multi_s.png b/deal.II/examples/step-24/doc/multi_s.png
new file mode 100644 (file)
index 0000000..9bbd989
Binary files /dev/null and b/deal.II/examples/step-24/doc/multi_s.png differ
diff --git a/deal.II/examples/step-24/doc/multi_s2.png b/deal.II/examples/step-24/doc/multi_s2.png
new file mode 100644 (file)
index 0000000..77ddcf5
Binary files /dev/null and b/deal.II/examples/step-24/doc/multi_s2.png differ
diff --git a/deal.II/examples/step-24/doc/multi_sf.png b/deal.II/examples/step-24/doc/multi_sf.png
new file mode 100644 (file)
index 0000000..55b0c09
Binary files /dev/null and b/deal.II/examples/step-24/doc/multi_sf.png differ
diff --git a/deal.II/examples/step-24/doc/one.png b/deal.II/examples/step-24/doc/one.png
new file mode 100644 (file)
index 0000000..f63dac9
Binary files /dev/null and b/deal.II/examples/step-24/doc/one.png differ
diff --git a/deal.II/examples/step-24/doc/one_movie.gif b/deal.II/examples/step-24/doc/one_movie.gif
new file mode 100644 (file)
index 0000000..edccbbf
Binary files /dev/null and b/deal.II/examples/step-24/doc/one_movie.gif differ
diff --git a/deal.II/examples/step-24/doc/one_s.png b/deal.II/examples/step-24/doc/one_s.png
new file mode 100644 (file)
index 0000000..934bf6a
Binary files /dev/null and b/deal.II/examples/step-24/doc/one_s.png differ
diff --git a/deal.II/examples/step-24/doc/one_s2.png b/deal.II/examples/step-24/doc/one_s2.png
new file mode 100644 (file)
index 0000000..784fb79
Binary files /dev/null and b/deal.II/examples/step-24/doc/one_s2.png differ
diff --git a/deal.II/examples/step-24/doc/one_sf.png b/deal.II/examples/step-24/doc/one_sf.png
new file mode 100644 (file)
index 0000000..f94482d
Binary files /dev/null and b/deal.II/examples/step-24/doc/one_sf.png differ
diff --git a/deal.II/examples/step-24/doc/project-1.tex b/deal.II/examples/step-24/doc/project-1.tex
new file mode 100644 (file)
index 0000000..23b3609
--- /dev/null
@@ -0,0 +1,91 @@
+\documentclass[english]{article}
+\usepackage[T1]{fontenc}
+\usepackage[latin1]{inputenc}
+\usepackage{geometry}
+\geometry{verbose,letterpaper,tmargin=1in,bmargin=1in,lmargin=1in,rmargin=1in}
+\usepackage{amsmath}
+
+\makeatletter
+\usepackage{babel}
+\makeatother
+\begin{document}
+
+\section{FEM model}
+
+\begin{quote}
+The problem is
+
+\begin{align*}
+\Delta\bar{p}-\frac{1}{c_{0}^{2}}\frac{\partial^{2}\bar{p}}{\partial^{2}t} & =f\end{align*}
+
+
+with initial condition 
+
+\begin{eqnarray*}
+\bar{p}(0,\mathbf{r}) & = & b(r)\end{eqnarray*}
+
+
+Let 
+
+\begin{eqnarray*}
+v & = & \frac{\partial\bar{p}}{\partial t}\end{eqnarray*}
+
+
+then we have
+
+\begin{alignat*}{1}
+\bar{p}_{t}-v & =0\\
+\Delta\bar{p}-\frac{1}{c_{0}^{2}}\, v_{t} & =f\end{alignat*}
+
+
+and absorbing boundary condition 
+
+\begin{eqnarray*}
+\frac{\partial\bar{p}}{\partial\mathbf{n}} & =- & \frac{1}{c_{0}}\frac{\partial\bar{p}}{\partial t}\end{eqnarray*}
+
+
+$\frac{\partial\bar{p}}{\partial\mathbf{n}}$ is the normal derivative
+at the boundary. This is a the time-varying FEM model. by discretizing
+according to $t$, we have
+
+\begin{eqnarray*}
+(\frac{\bar{p}^{n}-\bar{p}^{n-1}}{\delta t},\phi)_{\Omega}-\,(\theta v^{n}+(1-\theta)v^{n-1},\phi)_{\Omega} & = & 0\\
+-(\Delta((\theta\bar{p}^{n}+(1-\theta)\bar{p}^{n-1}),\bigtriangledown\phi)_{\Omega}-\frac{1}{c_{0}}(\frac{\bar{p}^{n}-\bar{p}^{n-1}}{\delta t},\phi)_{\partial\Omega}-\frac{1}{c_{0}^{2}}(\frac{v^{n}-v^{n-1}}{\delta t},\phi)_{\Omega} & = & (\theta f^{n}+(1-\theta)f^{n-1},\phi)_{\Omega}\end{eqnarray*}
+
+
+we obtain
+
+\begin{eqnarray*}
+M\bar{p}^{n}-(\delta t\,\theta)Mv^{n} & = & M\bar{p}^{n-1}+\delta t\,(1-\theta)\, M\, v^{n-1}\\
+(-c_{0}^{2}\,\delta t\,\theta A-c_{0}\, B)\bar{p}^{n}-Mv^{n} & = & (c_{0}^{2}\,\delta t\,(1-\theta)A-c_{0}B)\bar{p}^{n-1}-M\, v^{n-1}+c_{0}^{2}\delta t(\theta F^{n}+(1-\theta)F^{n-1})\end{eqnarray*}
+
+
+Write the above two equations as a matrix form
+
+\begin{eqnarray*}
+\left(\begin{array}{cc}
+M & -(\delta t\,\theta)M\\
+c_{0}^{2}\,\delta t\,\theta A+c_{0}\, B & M\end{array}\right)\left(\begin{array}{c}
+\bar{p}^{n}\\
+v^{n}\end{array}\right) & = & \left(\begin{array}{c}
+G_{1}\\
+G_{2}\end{array}\right)\end{eqnarray*}
+
+
+where 
+
+\begin{center}$\left(\begin{array}{c}
+G_{1}\\
+G_{2}\end{array}\right)=\left(\begin{array}{c}
+M\bar{p}^{n-1}+\delta t\,(1-\theta)Mv^{n-1}\\
+(-c_{0}^{2}\,\delta t\,(1-\theta)A+c_{0}B)\bar{p}^{n-1}+M\, v^{n-1}-c_{0}^{2}\delta t(\theta F^{n}+(1-\theta)F^{n-1})\end{array}\right)$\end{center}
+
+From the above matrix, we can obtain
+
+\begin{eqnarray*}
+(M+(\delta t\,\theta\, c_{0})^{2}A+c_{0}\,\delta t\,\theta\, B)\bar{p}^{n} & = & G_{1}+(\delta t\,\theta)G_{2}\\
+Mv^{n} & = & -(c_{0}^{2}\,\delta t\,\theta\, A+c_{0}B)\bar{p}^{n}+G_{2}\end{eqnarray*}
+
+\end{quote}
+
+\end{document}
diff --git a/deal.II/examples/step-24/doc/results.dox b/deal.II/examples/step-24/doc/results.dox
new file mode 100644 (file)
index 0000000..bc0a7f4
--- /dev/null
@@ -0,0 +1,90 @@
+<a name="Results"></a>
+<h1>Results</h1>
+
+The program writes values obtained by the detector to the disk. We then 
+draw them in plots. Experimental data are also collected for comparision. 
+Currently our experiments have only been done in two dimenstions by 
+circually scanning a single detector. We make the sample as a thin slice 
+in the X-Y plane (Z=0), and assume that signals from other Z directions 
+won't contribute to image construnction. By such an assumption we 
+only compare our experimental data with two dimensional simulated data.  
+
+
+<h3> One absorber </h3>
+
+This movie shows the thermoacoustic waves generated by a small absorber
+propagating in the medium (in our simulation, we assume the medium is mineral
+oil, which has a acoustic speed of 1.437 mm$\mu$s).
+
+@image html one_movie.gif
+
+The following two figures show the results obtained by the detector. The left
+side figure is obtained through experiments, and the right side figure is the 
+simulated data. In the experiment, a small strong absorber was embedded in weaker 
+absorbing tissue, we can still see the signals generated by the weak absoring
+tissue, which surrounds the signals generated by the small strong absorber in 
+the center. In the simulated data, we only simulate the small strong absorber.
+
+<TABLE WIDTH="100%">
+<tr>
+<td>
+@image html one.png
+</td>
+<td>
+@image html one_s.png
+</td>
+</tr>
+</table>
+
+In real situation, the detector has limited bandwidth. The thermoacoustic waves
+passing through the detector will be filtered. By using a high-pass filter,
+the simulated results look closer to the experimenal data. This is being done
+in MATLAB.
+
+@image html one_sf.png
+
+Because the mesh mismatch in the coarse mesh used in above the simulation, 
+we observe that waves are reflected before they travel to the boundary.
+The problem can be alleviated by using finer mesh. This result is not filtered.
+@image html one_s2.png
+
+
+<h3>Multiple absorbers</h3>
+
+We have done the simulation for one small absrober. To further verify the program,
+we will shown the simulation results for multiple absorbers.
+
+First we will use a movie to show the propagation of the generated thermoacoustic 
+waves in the medium by multiple absorbers:
+@image html multi_movie.gif
+
+Then the experimental data and our simulated data are compared in the following two
+figures:
+<TABLE WIDTH="100%">
+<tr>
+<td>
+@image html multi.png
+</td>
+<td>
+@image html multi_s.png
+</td>
+</tr>
+</table>
+The data points in experimental data are about eight times more than that in the 
+simulated data. In order to get similiar results as in the experiment, we need finer 
+mesh and smaller time step in the simulation. 
+Same as for the one absorber case, by using a high-pass filter, we get the following 
+results:
+
+@image html multi_sf.png
+
+By using finer mesh, we obtain better simulated results. The results shown is not filtered.
+
+@image html multi_s2.png
+
+
+
+
+
diff --git a/deal.II/examples/step-24/step-24 b/deal.II/examples/step-24/step-24
new file mode 100755 (executable)
index 0000000..061710c
Binary files /dev/null and b/deal.II/examples/step-24/step-24 differ
diff --git a/deal.II/examples/step-24/step-24.cc b/deal.II/examples/step-24/step-24.cc
new file mode 100644 (file)
index 0000000..e71de80
--- /dev/null
@@ -0,0 +1,697 @@
+
+/*    $Id: project.cc modified from heat-equation.cc 2006/03/05 $ */
+/*    Author: Xing Jin                                            */
+/*                                                                */
+/*    $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $        */
+/*    Version: $Name:  $                                          */
+/*                                                                */
+/*    Copyright (c) 1999,2000,2001,2002,2003,2004,2005,2006       */
+/*    by the deal.II authors.                                     */
+/*    This file is subject to QPL and may not be  distributed     */
+/*    without copyright and license information. Please refer     */
+/*    to the file deal.II/doc/license.html for the  text  and     */
+/*    further information on this license.                        */
+                                                          
+
+                           // @sect3{Include files}
+                           // Most include files have been covered in 
+                           // step-6 and will not be further commented on
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+                           // We will need to read the value at a specific 
+                           // location. This including file is needed for 
+                           // finding a cell that contains a given point
+#include <grid/grid_tools.h>
+#include <grid/grid_generator.h>
+                           // Because the scanning geometry is on a circle,
+                           // the boundaries are not straight lines, so
+                           // we need some classes to predefine some 
+                           // boundary description
+#include <grid/tria_boundary_lib.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+                         
+#include <fe/mapping_q1.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <dofs/dof_constraints.h>
+
+#include <numerics/matrices.h>
+#include <numerics/vectors.h>
+
+#include <numerics/data_out.h>
+                           // These are for c++
+#include <fstream>
+#include <iostream>
+#include <sstream>
+                                                                  
+#include <base/logstream.h>
+                             
+#include <base/point.h>
+
+                           // @sect3{"The forward problem" class template}
+
+                           // The main class is similar to the wave equation.
+                           // The difference is that we add an absorbing 
+                           // boundary condition. Because we are only interested
+                           // in values at specific locations, we define some
+                           // parameters to obtain the coordinates of those 
+                           // locations.
+template <int dim>
+class TATForwardProblem
+{
+  public:
+    TATForwardProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve_p ();
+    void solve_v ();
+    void output_results (const unsigned int timestep_number) const;
+
+    Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix_p;
+    SparseMatrix<double> system_matrix_v;
+    SparseMatrix<double> mass_matrix;
+    SparseMatrix<double> laplace_matrix;
+    SparseMatrix<double> boundary_matrix;
+                           // Number of refinement
+    unsigned int n_refinements;
+                           // The acoustic speed in the medium $c_0$
+    double acoustic_speed;
+                           // This parameter is needed for discritizing
+                           // time-dependent problem
+    double theta; 
+
+                           // The total data collection time
+    double total_time;
+                           // The size of the time step
+    double time_step;  
+                           // The detector circullarly scan the target region.
+                           // The step size of the detector is in angles
+    double step_angle;
+                           // The scanning radius
+    double radius;  
+   
+    
+    Vector<double>       solution_p;
+    Vector<double>       old_solution_p;
+    Vector<double>       system_rhs_p;
+
+    Vector<double>       solution_v;
+    Vector<double>       old_solution_v;
+    Vector<double>       system_rhs_v;        
+                                          
+      
+};
+
+                            // Declare a class template for the right hand side
+                            // of the pressure potential 
+template <int dim>
+class RightHandSide_p : public Function<dim> 
+{
+  public:
+    RightHandSide_p () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+                            // Declare a class template for the right hand side
+                            // of the derivative of the pressure potential                  
+template <int dim>
+class RightHandSide_v : public Function<dim> 
+{
+  public:
+    RightHandSide_v () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+                            // Declare a class template for the initial values
+                            // of the pressure potential
+template <int dim>
+class InitialValues_p : public Function<dim> 
+{
+  public:
+    InitialValues_p () : Function<dim>() {};
+    
+  virtual double value (const Point<dim> &p,
+                         const unsigned int  component = 0) const;
+};
+
+                            // Declare a class template for the initial values
+                            // of the derivative of the pressure potential
+template <int dim>
+class InitialValues_v : public Function<dim> 
+{
+  public:
+    InitialValues_v () : Function<dim>() {};
+    
+    virtual double value (const Point<dim> &p, 
+                         const unsigned int  component = 0) const;
+};
+
+                             // Here is the function to set the right hand side
+                             // values to be zero for pressure potential
+template <int dim>
+double RightHandSide_p<dim>::value (const Point<dim> &/*p*/,
+                                 const unsigned int /*component*/) const 
+{
+  return 0;
+}
+                              // Similarly we set the right-hand size of the 
+                              // derivative of the pressure potential to be 
+                              // zero
+template <int dim>
+double RightHandSide_v<dim>::value (const Point<dim> &/*p*/,
+                                 const unsigned int /*component*/) const 
+{
+  return 0;
+}
+
+
+                           // The sources of the thermoacoustic waves 
+                           // are small absorbers. We will compare the 
+                           // simulation results with the experimental
+                           // data.
+
+template <int dim>
+double InitialValues_p<dim>::value (const Point<dim> &p,
+                                  const unsigned int /*component*/) const       
+{
+                          
+  if (std::sqrt(p.square())< 0.025 )
+       return 1;
+                           // The "distance" function is used to compute
+                           // the Euclidian distance between two points.
+                                    
+  if (p.distance(Point<dim>(-0.135,0))<0.05)
+  return 1;
+                           
+  if (p.distance(Point<dim>(0.17,0))<0.03)
+   return 1;
+
+  if (p.distance(Point<dim>(-0.25,0))<0.02)
+      return 1;
+
+  if (p.distance(Point<dim>(-0.05,-0.15))<0.015)
+      return 1;
+
+    return 0;
+}
+                            // Initial value for the derivative of
+                            // pressure potential is set to zero
+template <int dim>   
+double InitialValues_v<dim>::value (const Point<dim> &/*p*/,
+                                  const unsigned int /*component*/) const 
+{
+ return 0;
+}
+
+                            // Evaluate point values at arbitrary locations 
+                            // In real situation, we collect data by placing
+                            // a detector in the medium. By scanning the detector,
+                            // we obtain a series projections of the target
+                            // from different viewing angles. By using a 
+                            // circular radon transform, we can reconstruct
+                            // the energy distribution in the target area from 
+                            // the measurements obtained by the detectors.
+                                 
+template <int dim> 
+double point_value (const DoFHandler<dim> &dof,
+            const Vector<double>  &fe_function,
+            const Point<dim>      &point)
+{                              
+                            // Define a map that maps the unit cell to a 
+                            // a general grid cell with straight lines in 
+                            // <dim> dimension
+  static const MappingQ1<dim> mapping;
+  const FiniteElement<dim>& fe = dof.get_fe();
+
+  Assert(fe.n_components() == 1,
+        ExcMessage ("Finite element is not scalar as is necessary for this function"));
+
+                             // First find the cell in which this point
+                             // is, initialize a quadrature rule with
+                             // it, and then a FEValues object
+                             // The algorithm will first look for the 
+                             // surrounding cell on a coarse grid, and
+                             // then recersively checking its sibling
+                             // cells.
+  const typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (dof, point);
+              
+  const Point<dim> unit_point = mapping.transform_real_to_unit_cell(cell, point);
+  Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
+          ExcInternalError());
+
+  const Quadrature<dim> quadrature (unit_point);
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell);
+
+                             // Then use this to get at the values of
+                             // the given fe_function at this point
+  std::vector<double> u_value(1);
+  fe_values.get_function_values(fe_function, u_value);
+
+  return u_value[0];
+}
+
+
+                             // @sect4{Initialize the problem}
+                             // Acoustic_speed here is the acoustic speed 
+                             // in the medium. Specifically we use acoustic speed
+                             // in mineral oil. We use Crank-Nicolson scheme
+                             // for our time-dependent problem, therefore theta is
+                             // set to be 0.5.  The step size of the detector
+                             // is 2.25 degree, which means we need 160 steps 
+                             // in order to finish a circular scan. The radius of the
+                             // scanning circle is select to be half way between  
+                             // the center and the boundary to avoid the reflections 
+                             // from the the boundary, so as to miminize the 
+                             // interference brought by the inperfect absorbing 
+                             // boundary condition. The time step is selected 
+                             // to satisfy $k = h/c$                           
+template <int dim>
+TATForwardProblem<dim>::TATForwardProblem () :
+                fe (1),
+               dof_handler (triangulation),
+               n_refinements (7),
+                acoustic_speed (1.437),
+                theta (0.5),
+               total_time (0.7),
+               time_step (0.5/std::pow(2.,1.0*n_refinements)/acoustic_speed),  
+                step_angle (2.25),  
+                radius (0.5)
+
+{}   
+
+                              // This is similar to step-6 except that
+                              // the mesh generated is a hyper_ball. We select
+                              // hyper_ball instead of hyper_cube because of
+                              // our data collection geometry is on a circular in 
+                              // 2-D, and on a sphere in 3-D. 
+template <int dim>
+void TATForwardProblem<dim>::make_grid_and_dofs ()
+{
+                              // In two dimensional domain. The center of the
+                              // circle shall be the point (0,0) and the radius
+                              // is 1
+  const Point<2> center (0,0);
+  GridGenerator::hyper_ball (triangulation, center, 1);
+                              // Accordingly, we use hyper ball boundary
+                              // instead of hyper cube.
+  static const HyperBallBoundary<dim> boundary_description(center);
+  triangulation.set_boundary (0,boundary_description);
+                              //  The mesh is refined n_refinements times
+  triangulation.refine_global (n_refinements);
+
+  std::cout << "   Number of active cells: "
+           << triangulation.n_active_cells()
+           << std::endl
+           << "   Total number of cells: "
+           << triangulation.n_cells()
+           << std::endl;
+
+  dof_handler.distribute_dofs (fe);
+
+  std::cout << "   Number of degrees of freedom: "
+           << dof_handler.n_dofs()
+           << std::endl;
+
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+                              // We will do the following for both
+                              // the pressure potential and its derivative
+  system_matrix_p.reinit (sparsity_pattern);
+  system_matrix_v.reinit (sparsity_pattern);
+
+  mass_matrix.reinit (sparsity_pattern);
+  laplace_matrix.reinit (sparsity_pattern);
+  boundary_matrix.reinit (sparsity_pattern);
+
+  solution_p.reinit (dof_handler.n_dofs());
+  old_solution_p.reinit (dof_handler.n_dofs());
+  system_rhs_p.reinit (dof_handler.n_dofs());
+
+  solution_v.reinit (dof_handler.n_dofs());
+  old_solution_v.reinit (dof_handler.n_dofs());
+  system_rhs_v.reinit (dof_handler.n_dofs());
+}
+
+
+                               // @sect3{ Assemble system}
+                               // Because we used absorbing boundary condition in the
+                               // simulation, a new boundary matrix is introduced.
+                               // We need to assemble boundary matrix. The detailed
+                               // description for assembling matrix is discussed in
+                               // step-3. 
+template <int dim>
+void TATForwardProblem<dim>::assemble_system () 
+{  
+  MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3),
+                                    mass_matrix);
+  MatrixCreator::create_mass_matrix (dof_handler, QTrapez<dim>(),
+                                    mass_matrix);
+  mass_matrix /= 2;
+
+  MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3),
+                                       laplace_matrix);
+  MatrixCreator::create_laplace_matrix (dof_handler, QTrapez<dim>(),
+                                       laplace_matrix);
+  laplace_matrix /= 2;
+
+  const QGauss<dim-1>  quadrature_formula(3);
+
+  FEFaceValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values  |  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);
+
+  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)
+    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+      if (cell->at_boundary(f))
+       {
+         cell_matrix = 0;
+
+         fe_values.reinit (cell, f);
+
+         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) += (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)
+             boundary_matrix.add (local_dof_indices[i],
+                                local_dof_indices[j],
+                                cell_matrix(i,j));
+       }
+
+                                 // The system matrix of pressure potential 
+                                 // is shown in the introduction 
+  system_matrix_p = 0;
+  system_matrix_p.copy_from (mass_matrix);
+  system_matrix_p.add (time_step*time_step*theta*theta*acoustic_speed*acoustic_speed, laplace_matrix);
+  system_matrix_p.add (acoustic_speed*theta*time_step, boundary_matrix);
+                                 // The system matrix of the derivative
+                                 // of the pressure potential is same as 
+                                 // the mass matrix
+  system_matrix_v = 0;
+  system_matrix_v.copy_from (mass_matrix);
+
+}
+
+
+                                 // We will solve two equations. 
+                                 // We first solve for pressure potential
+                                 // at a time step, then solve the derivative
+                                 // of the pressure potential.
+template <int dim>
+void TATForwardProblem<dim>::solve_p () 
+{
+  SolverControl           solver_control (1000, 1e-10);
+  SolverCG<>              cg (solver_control);
+  cg.solve (system_matrix_p, solution_p, system_rhs_p,
+           PreconditionIdentity());
+
+  std::cout << "   " << solver_control.last_step()
+           << " CG iterations needed to obtain convergence."
+           << std::endl;
+}
+                       
+                                 // To solve the derivative of the pressure potential
+template <int dim>
+void TATForwardProblem<dim>::solve_v () 
+{
+  SolverControl           solver_control (1000, 1e-10);
+  SolverCG<>              cg (solver_control);
+  cg.solve (system_matrix_v, solution_v, system_rhs_v,
+           PreconditionIdentity());
+
+  std::cout << "   " << solver_control.last_step()
+           << " CG iterations needed to obtain convergence."
+           << std::endl;
+}
+
+                                 // We output the solution for pressure potential
+                                 // at each time step in "vtk" format.
+template <int dim>
+void TATForwardProblem<dim>::output_results (const unsigned int timestep_number) const
+{
+  
+  DataOut<dim> data_out; 
+  
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution_p, "P");
+  data_out.add_data_vector (solution_v, "V");
+
+  data_out.build_patches ();
+
+  std::ostringstream filename;
+  filename << "solution-"
+          << timestep_number<<".vtk";
+  std::ofstream output (filename.str().c_str());
+  data_out.write_vtk (output);
+
+  
+}
+
+                                 // This is the main function
+                                 // 
+                                      
+template <int dim>
+void TATForwardProblem<dim>::run () 
+{
+  std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
+  
+  make_grid_and_dofs();
+  assemble_system ();
+
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  VectorTools::project (dof_handler,constraints, 
+                            QGauss<dim>(3), InitialValues_p<dim>(),
+                            old_solution_p);
+  VectorTools::project (dof_handler,constraints, 
+                            QGauss<dim>(3), InitialValues_v<dim>(),
+                            old_solution_v);
+
+
+  unsigned int timestep_number = 1;
+  unsigned int n_steps;
+  unsigned int n_detectors;
+  double scanning_angle;
+
+                                // Number of time steps is defined as the
+                                // ratio of the total time to the time step                                 
+  n_steps=static_cast<unsigned int>(std::floor(total_time/time_step));     
+                                // Number of detector positions is defined          
+                                // as the ratio of 360 degrees to the step
+                                // angle
+  n_detectors=static_cast<unsigned int>(std::ceil(360/step_angle));
+                                // Define two vectors to hold the coordinates
+                                // of the detectors in the scanning
+                                // geometry
+  Vector<double> detector_x (n_detectors+1);
+  Vector<double> detector_y (n_detectors+1);
+                                // Define a vector to hold the value obtained
+                                // by the detector
+  Vector<double> project_dat (n_steps * n_detectors +1);
+                                // Get the coordinates of the detector on the 
+                                // different locations of the circle.
+                                // Scanning angle is viewing angle at 
+                                // current position. The coordinates of
+                                // the detectors are calculated from the radius
+                                // and scanning angle.
+  scanning_angle=0;
+  for (unsigned int i=1; i<=n_detectors;  i++){
+                                // Scanning clockwisely. We need to change the angles
+                                // into radians because std::cos and std:sin accept
+                                // values in radian only
+    scanning_angle -= step_angle/180 * 3.14159265;   
+    detector_x(i) = radius * std::cos(scanning_angle);
+    detector_y(i) = radius * std::sin(scanning_angle);
+  }
+
+  std::cout<< "Total number of time steps = "<< n_steps <<std::endl;
+  std::cout<< "Total number of detectors = "<< n_detectors << std::endl;
+
+                               // Open a file to write the data
+                               // obtained by the detectors 
+                         
+  std::ofstream proj_out;
+  proj_out.open("proj.dat");
+  
+
+  for (double time = time_step; time<=total_time; time+=time_step, ++timestep_number)
+    {
+      std::cout << std::endl;                                       
+      std::cout<< "time_step " << timestep_number << " @ t=" << time << std::endl;
+
+      Vector<double> tmp1 (solution_p.size());
+      Vector<double> tmp2 (solution_v.size());
+      Vector<double> F1 (solution_p.size());
+      Vector<double> F2 (solution_v.size());
+                               // Calculate G1 as defined in the introduction section
+                         
+      mass_matrix.vmult (tmp1, old_solution_p); 
+      mass_matrix.vmult (tmp2, old_solution_v); 
+      F1 = tmp1;
+      F1.add(time_step * (1-theta), tmp2); 
+                               // Calculate G2 as defined in the introduction section
+      mass_matrix.vmult (tmp1, old_solution_v);
+      laplace_matrix.vmult (tmp2, old_solution_p); 
+      F2 = tmp1;
+      F2.add(-acoustic_speed*acoustic_speed*time_step*(1-theta), tmp2);
+      tmp1=0;
+      boundary_matrix.vmult (tmp1,old_solution_p);
+      F2.add(acoustic_speed,tmp1);
+      
+                               // Compute the pressure potential p, the formula
+                               // has been presented in the introduction section
+
+      system_rhs_p = F1; 
+      system_rhs_p.add(time_step * theta , F2);
+
+      RightHandSide_p<dim> rhs_function_p;
+      rhs_function_p.set_time (time);
+
+      tmp1=0;
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                           rhs_function_p, tmp1);
+    
+      system_rhs_p.add(-theta * theta * time_step * time_step*acoustic_speed*acoustic_speed,tmp1); 
+      rhs_function_p.set_time (time-time_step);
+      tmp1=0;
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                           rhs_function_p, tmp1);
+      
+      system_rhs_p.add(-theta * (1-theta) * time_step * time_step*acoustic_speed*acoustic_speed,tmp1); 
+
+      solve_p ();
+
+                                // Compute the derivative potential pressure.
+                                // The formula has been presented in the introduction
+                                // section. The potential derivative is calculated 
+                                // after the potential pressure because the calculation
+                                // depends on the current value of the potential 
+                                // pressure
+
+      system_rhs_v = F2;
+      tmp1 = 0;
+      laplace_matrix.vmult (tmp1, solution_p);
+      system_rhs_v.add(-time_step * theta*acoustic_speed*acoustic_speed, tmp1);
+      tmp1 = 0;
+      boundary_matrix.vmult(tmp1, solution_p);
+      system_rhs_v.add(-acoustic_speed,tmp1);
+      
+      RightHandSide_v<dim> rhs_function_v;
+      rhs_function_v.set_time (time); 
+
+      tmp2 = 0;
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                           rhs_function_v, tmp2);
+
+      system_rhs_p.add(-theta * time_step*acoustic_speed*acoustic_speed,tmp2); 
+
+      rhs_function_v.set_time (time-time_step);
+      tmp2 = 0;
+      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
+                                           rhs_function_v, tmp2);
+      system_rhs_p.add(-(1-theta)*time_step*acoustic_speed*acoustic_speed,tmp2);
+      
+      solve_v ();
+                               // Compute the energy in the system.By checking
+                               // energy change in the system, we can verify
+                               // the correctness of the code. 
+
+      double energy = (mass_matrix.matrix_scalar_product(solution_v,solution_v)+
+                   acoustic_speed*acoustic_speed*laplace_matrix.matrix_scalar_product(solution_p,solution_p))/2;        
+                                                                         
+      std::cout << "energy= " << energy << std::endl;
+
+      //  output_results (timestep_number);
+      
+                               //  Evaluate the value at specific locations. 
+                               //  For 2-D, it is on a circle. For 1-D, 
+                               //  it is a point detector.
+
+      proj_out << time ;
+
+      for (unsigned i=1 ; i<=n_detectors; i++){
+        project_dat((timestep_number-1)*n_detectors+i)=point_value (dof_handler,solution_p, 
+                                                   Point<2>(detector_x(i),detector_y(i)));
+        proj_out << " "<< project_dat((timestep_number-1)*n_detectors+i)<<" " ;
+      }
+
+      proj_out<<std::endl;
+          
+                              // Update the values for the pressure potential 
+                              // and its derivative. 
+         
+      old_solution_p = solution_p;
+      solution_p = 0;
+      old_solution_v = solution_v;      
+      solution_v = 0;
+    }
+      proj_out.close();
+                             
+    
+}
+                              // @sect3{The "main" function}
+                              // The main function calls the above functions
+                              // in the order of their appearances.
+
+int main ()  
+{
+  deallog.depth_console (0);
+  {
+    TATForwardProblem<2> TAT_forward_2d;
+    TAT_forward_2d.run ();
+  }
+  
+  return 0;
+}
diff --git a/deal.II/examples/step-25/Makefile b/deal.II/examples/step-25/Makefile
new file mode 100644 (file)
index 0000000..f291d79
--- /dev/null
@@ -0,0 +1,154 @@
+# $Id: Makefile 11909 2005-12-21 13:30:36Z guido $
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+################################################################
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<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-25/Makefile.dep b/deal.II/examples/step-25/Makefile.dep
new file mode 100644 (file)
index 0000000..4ac559b
--- /dev/null
@@ -0,0 +1,244 @@
+./step-25.o: \
+               step-25.cc\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/config.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/data_out_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function_time.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/geometry_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/logstream.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/multithread_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/numbers.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/point.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/polynomial.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/smartpointer.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/subscriptor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_product_polynomials.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/thread_management.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/vector_slice.h\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/assert.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/checked_delete.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/posix_features.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_compiler_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_platform_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_stdlib_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/suffix.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/current_function.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_gcc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_win32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/bad_weak_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/interlocked.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lightweight_mutex.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_nop.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_win32_cs.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/quick_allocator.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_ptr_nmt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_cw_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ia64.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_x86.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_nt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_pt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_w32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_impl.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/workaround.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/shared_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/throw_exception.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/tuple/tuple.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/alignment_of.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/type_with_alignment.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_constraints.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/function_map.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_collection.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_poly.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_q.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_update_flags.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_values.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_generator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_faces.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_hex.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_line.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_quad.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out_stack.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/matrices.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/vectors.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/full_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/identity_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/precondition.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_control.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_gmres.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparse_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparsity_pattern.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector_memory.h
+./step-25.g.o: \
+               step-25.cc\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/config.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/data_out_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/function_time.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/geometry_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/logstream.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/multithread_info.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/numbers.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/point.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/polynomial.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/quadrature_lib.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/smartpointer.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/subscriptor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/table_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/tensor_product_polynomials.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/thread_management.h\
+               /home/bangerth/p/deal.II/1/deal.II/base/include/base/vector_slice.h\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/assert.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/checked_delete.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/posix_features.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_compiler_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_platform_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/select_stdlib_config.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/config/suffix.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/current_function.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_gcc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/atomic_count_win32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/bad_weak_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/interlocked.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lightweight_mutex.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_nop.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_pthreads.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/lwm_win32_cs.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/quick_allocator.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_count.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/shared_ptr_nmt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_cw_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ia64.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_ppc.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_gcc_x86.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_nt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_pt.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_base_w32.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/sp_counted_impl.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/detail/workaround.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/shared_ptr.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/throw_exception.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/tuple/tuple.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/alignment_of.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/contrib/boost/include/boost/type_traits/type_with_alignment.hpp\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_constraints.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/dof_tools.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/function_map.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/dofs/hp_dof_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_collection.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_poly.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_q.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_update_flags.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/fe_values.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/fe/mapping.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/grid_generator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_accessor.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_faces.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_hex.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator.templates.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_levels.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_line.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_objects.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/grid/tria_quad.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_accessor.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_handler.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/multigrid/mg_dof_iterator_selector.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/data_out_stack.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/matrices.h\
+               /home/bangerth/p/deal.II/1/deal.II/deal.II/include/numerics/vectors.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_indices.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/block_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/exceptions.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/full_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/identity_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_block_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_parallel_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/petsc_vector_base.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/precondition.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_control.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/solver_gmres.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparse_matrix.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/sparsity_pattern.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector.h\
+               /home/bangerth/p/deal.II/1/deal.II/lac/include/lac/vector_memory.h
diff --git a/deal.II/examples/step-25/doc/animate.sh b/deal.II/examples/step-25/doc/animate.sh
new file mode 100755 (executable)
index 0000000..31807c2
--- /dev/null
@@ -0,0 +1,38 @@
+#!/bin/sh
+dim=1
+tstep=0
+tstepinc=20
+numtsteps=6300
+batchfile='animation.plt'
+
+# optional clean up...all these files will get overwritten
+rm "solution-"$dim"d-"*".png"
+rm "solution-"$dim"d.gif"
+
+# always gotta delete this one, though
+rm $batchfile
+
+# generate the gnuplot batch script to plot all the desire time steps
+while [ $tstep -lt $numtsteps ]
+do
+    if [ $tstep -lt 10 ]; then
+       ztstep="0000"$tstep
+    elif [ $tstep -lt 100 ]; then
+       ztstep="000"$tstep
+    elif [ $tstep -lt 1000 ]; then
+       ztstep="00"$tstep
+    elif [ $tstep -lt 10000 ]; then
+       ztstep="0"$tstep
+    else
+       ztstep=$tstep
+    fi
+
+    echo "call \"plot.plt\" $dim $ztstep" >> $batchfile
+    
+    let tstep=tstep+$tstepinc
+done
+
+gnuplot -persist $batchfile
+
+# use ImageMagick to create an animated gif from the PNG files
+convert -delay 0 -loop 0 "solution-"$dim"d-*.png" "solution-"$dim"d.gif"
\ No newline at end of file
diff --git a/deal.II/examples/step-25/doc/christov.1d-breather.gif b/deal.II/examples/step-25/doc/christov.1d-breather.gif
new file mode 100644 (file)
index 0000000..9d38a39
Binary files /dev/null and b/deal.II/examples/step-25/doc/christov.1d-breather.gif differ
diff --git a/deal.II/examples/step-25/doc/christov.1d-breather_stp.png b/deal.II/examples/step-25/doc/christov.1d-breather_stp.png
new file mode 100644 (file)
index 0000000..09a8728
Binary files /dev/null and b/deal.II/examples/step-25/doc/christov.1d-breather_stp.png differ
diff --git a/deal.II/examples/step-25/doc/christov.2d-angled_kink.gif b/deal.II/examples/step-25/doc/christov.2d-angled_kink.gif
new file mode 100644 (file)
index 0000000..5fe6b63
Binary files /dev/null and b/deal.II/examples/step-25/doc/christov.2d-angled_kink.gif differ
diff --git a/deal.II/examples/step-25/doc/christov.2d-kink.gif b/deal.II/examples/step-25/doc/christov.2d-kink.gif
new file mode 100644 (file)
index 0000000..d2df430
Binary files /dev/null and b/deal.II/examples/step-25/doc/christov.2d-kink.gif differ
diff --git a/deal.II/examples/step-25/doc/christov.2d-kink.png b/deal.II/examples/step-25/doc/christov.2d-kink.png
new file mode 100644 (file)
index 0000000..2782629
Binary files /dev/null and b/deal.II/examples/step-25/doc/christov.2d-kink.png differ
diff --git a/deal.II/examples/step-25/doc/christov.2d-pseudobreather.gif b/deal.II/examples/step-25/doc/christov.2d-pseudobreather.gif
new file mode 100644 (file)
index 0000000..1c6826d
Binary files /dev/null and b/deal.II/examples/step-25/doc/christov.2d-pseudobreather.gif differ
diff --git a/deal.II/examples/step-25/doc/intro.dox b/deal.II/examples/step-25/doc/intro.dox
new file mode 100644 (file)
index 0000000..0aa5ab5
--- /dev/null
@@ -0,0 +1,87 @@
+<a name="Intro"></a> <h1>Introduction</h1>
+The goal of this program is to solve the sine-Gordon soliton equation in 1, 2 or 3 spatial dimensions. The motivation for solving this equation is that very little is known about the nature of the solutions in 2D and 3D, even though the 1D case has been studied extensively.
+
+Rather facetiously, the sine-Gordon equation's moniker is a pun on the so-called Klein-Gordon equation, which is the relativistic version of the Schr&ouml;dinger equation. The resemblence is not just superficial, the sine-Gordon equation has been shown to model some unified-field phenomena such as interaction of subatomic particles (see, e.g., Perring &amp; Skyrme in Nuclear Physics <b>31</b>) and the Josephson (quantum) effect in semiconductors junctions (see, e.g., <a href="http://en.wikipedia.org/wiki/Long_Josephson_junction">http://en.wikipedia.org/wiki/Long_Josephson_junction</a>). Furthermore, from the mathematical standpoint, since the sine-Gordon equation is "completely integrable," it is a candidate for study using the usual methods such as the inverse scattering transform. Consequently, over the years, many interesting solitary-wave, and even stationary, solutions to the sine-Gordon equation have been found. For more on the sine-Gordon equation, the inverse scattering transform and other methods for finding analyical soliton equations, the reader should consult the following "classical" references on the subject: G. L. Lamb's <i>Elements of Soliton Theory</i> (Chapter 5, Section 2) and G. B. Whitham's <i>Linear and Nonlinear Waves</i> (Chapter 17, Sections 10-13).
+
+<h3>Statement of the problem</h3>
+The sine-Gordon initial-boundary-value problem (IBVP) we wish to solve is
+\f{eqnarray*}
+  u_{tt}-\Delta u &=& -\sin(u) \quad\mbox{for}\quad (x,t) \in \Omega \times (t_0,t_f],\\
+  {\mathbf n} \cdot \nabla u &=& 0 \quad\mbox{for}\quad (x,t) \in \partial\Omega 
+          \times (t_0,t_f],\\
+  u(x,t_0) &=& u_0(x).
+\f}
+We have chosen to enforce zero Neumann boundary conditions in order for waves to reflect off the boundaries of our domain. It should be noted, however, that Dirichlet boundary conditions are not appropriate for this problem. Even though the solutions to the sine-Gordon equation are localized, it only makes sense to specify (Dirichlet) boundary conditions at $x=\pm\infty$, otherwise either a solution does not exist or only the trivial solution $u=0$ exists.
+
+However, the form of the equation above is not ideal for numerical discretization. If we were to discretize the second-order time derivative directly and accurately, then  we would need a large stencil (i.e., several time steps would need to kept in the memory), which could become expensive. Therefore, we split the second-order (in time) sine-Gordon equation into a system of two first-order (in time) equations, which we call the split, or velocity, formulation. To this end, by setting $v = u_t$, it is easy to see that the sine-Gordon equation is equivalent to
+\f{eqnarray*}
+  u_t - v &=& 0,\\
+  v_t - \Delta u &=& -\sin(u). 
+\f}
+
+<h3>Discretization of the equations in time</h3>
+Now, we can discretize the split formulation in time using the the $\theta$-method, which has a stencil of only two time steps. By choosing a $\theta\in [0,1]$, the latter discretization allows us to choose from a continuum of schemes. In particular, if we pick $\theta=0$ or $\theta=1$, we obtain the first-order accurate explicit or implicit Euler method, respectively. Another important choice is $\theta=\frac{1}{2}$, which gives the second-order accurate Crank-Nicolson scheme. Henceforth, a superscript $n$ denotes the values of the variables at the $n^{\mathrm{th}}$ time step, i.e. at $t=t_n\equiv n k$, where $k$ is (fixed) the time step size. Thus, the split formulation of the sine-Gordon equation becomes
+\f{eqnarray*}
+  \frac{u^n - u^{n-1}}{k} - \left[\theta v^n + (1-\theta) v^{n-1}\right] &=& 0,\\
+  \frac{v^n - v^{n-1}}{k} - \Delta\left[\theta u^n + (1-\theta) u^{n-1}\right]
+  &=& -\sin\left[\theta u^n + (1-\theta) u^{n-1}\right]. 
+\f}
+
+We can simplify the latter via a bit of algebra. Eliminating $v^n$ from the first equation and rearranging, we obtain
+\f{eqnarray*}
+  \left[ 1-k^2\theta^2\Delta \right] u^n &=&
+        \left[ 1+k^2\theta(1-\theta)\Delta\right] u^{n-1} + k v^{n-1} 
+        - k^2\theta\sin\left[\theta u^n + (1-\theta) u^{n-1}\right],\\
+   v^n &=& v^{n-1} + k\Delta\left[ \theta u^n + (1-\theta) u^{n-1}\right]
+                - k\sin\left[ \theta u^n + (1-\theta) u^{n-1} \right].
+\f}
+
+It may seem as though we can just proceed to discretize the equations in space at this point, however, that would not work for all $\theta$ since the first equation above is nonlinear. Therefore, a nonlinear solver must be implemented, then equations can be discretized in space and solved. 
+
+To this end, we can use Newton's method. Given the nonlinear equation $F(u^n) = 0$, we produce successive approximations to $u^n$ as follows:
+\f{eqnarray*}
+  \mbox{ Find } \delta u^n_l \mbox{ s.t. } F'(u^n_l)\delta u^n_l = -F(u^n_l) 
+  \mbox{, set }  u^n_{l+1} = u^n_l + \delta u^n_l.
+\f}
+The iteration can be initialized with the old time step, i.e. $u^n_0 = u^{n-1}$, and eventually it will produce a solution to the first equation of the split formulation (see above). For the time discretizaion of the sine-Gordon under consideration here, we have that 
+\f{eqnarray*}
+  F(u^n_l) &=&  \left[ 1-k^2\theta^2\Delta \right] u^n_l -  
+                \left[ 1+k^2\theta(1-\theta)\Delta\right] u^{n-1} - k v^{n-1} 
+                + k^2\theta\sin\left[\theta u^n_l + (1-\theta) u^{n-1}\right],\\
+  F'(u^n_l) &=& 1-k^2\theta^2\Delta - k^2\theta^2\cos\left[\theta u^n_l 
+                       + (1-\theta) u^{n-1}\right].
+\f}
+Notice that while $F(u^n_l)$ is a function, $F'(u^n_l)$ is an operator.
+
+<h3>Weak formulation of the time-discretized equations</h3>
+With hindsight, we choose both the solution and the test space to be $H^1(\Omega)$. Hence, multiplying by a test function $\varphi$ and integrating, we obtain the following variational (or weak) formulation of the split formulation (including the nonlinear solver for the first equation) at each time step:
+\f{eqnarray*}
+  &\mbox{ Find}& \delta u^n_l \in H^1(\Omega) \mbox{ s.t. } 
+  \left( F'(u^n_l)\delta u^n_l, \varphi \right)_{\Omega} 
+  = -\left(F(u^n_l), \varphi \right)_{\Omega} \;\forall\varphi\in H^1(\Omega),
+  \mbox{ set } u^n_{l+1} = u^n_l + \delta u^n_l,\; u^n_0 = u^{n-1}.\\
+  &\mbox{ Find}& v^n \in H^1(\Omega) \mbox{ s.t. }
+  \left( v^n, \varphi \right)_{\Omega} = \left( v^{n-1}, \varphi \right)_{\Omega} 
+        - k\theta\left( \nabla u^n, \nabla\varphi \right)_{\Omega} 
+        - k (1-\theta)\left( \nabla u^{n-1}, \nabla\varphi \right)_{\Omega} 
+        - k\left(\sin\left[ \theta u^n + (1-\theta) u^{n-1} \right],
+        \varphi \right)_{\Omega} \;\forall\varphi\in H^1(\Omega).
+\f}
+Note that the we have used integration by parts and the zero Neumann boundary conditions on all terms involving the Laplacian operator ($\Delta$). Moreover, $F(\cdot)$ and $F'(\cdot)$ are as defined above, and $(\cdot,\cdot)_{\Omega}$ denotes the usual $L^2$ inner product over the domain $\Omega$, i.e. $(f,g)_{\Omega} = \int_\Omega fg \,\mathrm{d}x$. Finally, notice that the first equation is, in fact, the definition of an interative procedure, so it solved multiple times until a stopping criterion is met.
+
+<h3>Discretization of the weak formulation in space</h3>
+Using the Finite Element Method, we discretize the variation formulation in space. To this end, let $V_h$ be a finite-dimensional $H^1(\Omega)$-conforming finite element space ($\mathrm{dim}\, V_h = N < \infty$) with nodal basis $\{\varphi_1,\ldots,\varphi_N\}$. Hence, we can expand all functions in the weak formulation (see above) in terms of the nodal basis. Henceforth, we shall denote by a capital letter the vector of coefficients (in the nodal basis) of a function denoted by the same letter in lower case; e.g., $u^n = \sum_{i=1}^N U^n_i \varphi_i$ where $U^n \in {R}^N$ and $u^n \in H^1(\Omega)$. Thus, the finite-dimensional version of variation fomulation requires that we solve the following matrix equations at each time step:
+@f{eqnarray*}
+  F_h'(U^n_l)\delta U^n_l &=& -F_h(U^n_l), \quad
+       U^n_{l+1} = U^n_l + \delta U^n_l, \quad U^n_0 = U^{n-1}; \\
+  MV^n &=& MV^{n-1} - k \theta AU^n -k (1-\theta) AU^{n-1} - k S(u^n,u^{n-1}).
+@f}
+Above, the matrix $F_h'(\cdot)$ and the vector $F_h(\cdot)$ denote the discrete versions of the gadgets discussed above, i.e.
+\f{eqnarray*}
+  F_h(U^n_l) &=&  \left[ M+k^2\theta^2A \right] U^n_l -  
+               \left[ M-k^2\theta(1-\theta)A \right] U^{n-1} - k MV^{n-1} 
+               + k^2\theta S(u^n_l, u^{n-1}),\\
+  F_h'(U^n_l) &=& M+k^2\theta^2A 
+                               - k^2\theta^2N(u^n_l,u^{n-1})
+\f}
+Again, note that the first matrix equation above is, in fact, the defition of an iterative procedure, so it is solve multiple times until a stopping criterion is met. Moreover, $M$ is the mass matrix, i.e. $M_{ij} = \left( \varphi_i,\varphi_j \right)_{\Omega}$, $A$ is the Laplace matrix, i.e. $A_{ij} = \left( \nabla \varphi_i, \nabla \varphi_j \right)_{\Omega}$, $S$ is the nonlinear term in the auxilliary equation, i.e. $S_j(f,g) = \left( \sin\left[ \theta f + (1-\theta) g\right], \varphi_j \right)_{\Omega}$, and $N$ is the nonlinear term in the Jacobian matrix of $F(\cdot)$, i.e. $N_{ij}(f,g) = \left( \cos\left[ \theta f + (1-\theta) g\right]\varphi_i, \varphi_j \right)_{\Omega}$.
diff --git a/deal.II/examples/step-25/doc/plot.plt b/deal.II/examples/step-25/doc/plot.plt
new file mode 100644 (file)
index 0000000..2c63689
--- /dev/null
@@ -0,0 +1,20 @@
+set data style lines
+set hidden3d
+
+set xrange [-10:10]
+set yrange [-2.2:2.2]
+set zrange [-1:8]
+set xlabel "x"
+if ($0==1) set ylabel "u"
+if ($0==2) set ylabel "y"
+if ($0==2) set zlabel "u"
+
+set terminal png
+set output "solution-$0d-$1.png" 
+
+if ($0==1) pl "solution-$0d-$1.gpl"
+if ($0==2) spl "solution-$0d-$1.gpl"
+
+#set terminal x11
+#replot
+#pause 5
diff --git a/deal.II/examples/step-25/doc/results.dox b/deal.II/examples/step-25/doc/results.dox
new file mode 100644 (file)
index 0000000..c7d9951
--- /dev/null
@@ -0,0 +1,47 @@
+<a name="Results"></a> <h1>Results</h1>
+The explicit Euler time stepping scheme  ($\theta=0$) performs adequately for the problems we wish to solve. Unfortunately, a rather small time step has to be chosen due to stability issues --- $k\sim h/10$ appears to work for most the simulations we performed. On the other hand, the Crank-Nicolson scheme ($\theta=\frac{1}{2}$) is unconditionally stable, and (at least for the case of the 1D breather) we can pick the time step to be as large as $25h$ without any ill effects on the solution. The impicit Euler scheme ($\theta=1$) is "exponentially damped," so it is not a good choice for solving the sine-Gordon equation, which is conservative. However, some of the damped schemes in the continuum that is offered by the $\theta$-method were useful for eliminating spurious oscillations due to boundary effects.
+
+In the simulations below, we solve the sine-Gordon on the interval $\Omega = [-10,10]$ in 1D and on the square $\Omega = [-10,10]\times [-10,10]$ in 2D. In each case, the respective grid is refined uniformly 6 times ($\Rightarrow h\sim 2^{-6}$).
+
+<h3>An (1+1)D Solution</h3>
+The first example we discuss is the so-called 1D (stationary) breather solution of the sine-Gordon equation. The breather has the following closed-form expression:
+\f[
+u_{\mathrm{breather}}(x,t) = -4\arctan \left(\frac{m}{\sqrt{1-m^2}} \frac{\sin\left(\sqrt{1-m^2}t +c_2\right)}{\cosh(mx+c_1)} \right), 
+\f]
+where $c_1$, $c_2$ and $m<1$ are constants. In the simulation below, we have chosen $c_1=0$, $c_2=0$, $m=0.5$. Moreover, it is know that the period of oscillation of the breather is $2\pi\sqrt{1-m^2}$, hence we have chosen $t_0=-5.4414$ and $t_f=2.7207$ so that we can observe three oscillations of the solution. Then, taking $u_0(x) = u_{\mathrm{breather}}(x,t_0)$, $\theta=0$ and $k=h/10$, the program computed the following solution.
+
+\image html christov.1d-breather.gif "Animation of the 1D stationary breather." width=5cm
+
+Another way to visualize the (1+1)D solution is to use the output from the DataOutStack object to produce the space-time plot below instead of the animation above.
+
+\image html christov.1d-breather_stp.png "A space-time plot of the 1D stationary breather." width=5cm
+
+Furthermore, since the breather is an analytical solution of the sine-Gordon equation, we can use it to validate our code, at least in 1D. For the simulation shown in the two images above, the $L^2$ norm of the error in the finite element solution at each time step remained on the order of $10^{-2}$. Hence, we can conclude that the numerical method has been implemented correctly in the program.
+
+<h3>A few (2+1)D Solutions</h3>
+The only analytical solution to the sine-Gordon equation in (2+1)D that can be found in the literature is the so-called kink solitary wave. It has the following closed-form expression:
+\f[
+u_{\mathrm{kink}}(x,t) = -4\arctan \left\{a_0\left[x\cos\vartheta+\sin\vartheta\left(y\cosh\lambda + t\sinh\lambda\right) \right]\right\},
+\f]
+where $a_0$, $\vartheta$ and $\lambda$ are constants. In the simulation below we have chosen $a_0=\lambda=1$. Notice that if $\vartheta=\pi$ the kink is stationary, hence it would make a good solution against which we can validate the program in 2D because no reflections off the boundary of the domain occur.
+
+The simulation shown below was performed with $u_0(x) = u_{\mathrm{kink}}(x,t_0)$, $\theta=\frac{1}{2}$, $k=20h$, $t_0=1$ and $t_f=500$. The $L^2$ norm of the error of the finite element solution at each time step remained on the order of $10^{-2}$, showing that the program is working correctly in 2D, as well as 1D. Unfortunately, the solution is not very interesting, nonetheless we have included a snapshot of it below for completeness. 
+
+\image html christov.2d-kink.png "Stationary 2D kink." width=5cm
+
+Now that we have validated the code in 1D and 2D, we move to a problem where an analytical solution is unknown. 
+
+To this end, we can rotate the kink solution discussed above about the $z$ axis, e.g. let  $\vartheta=\frac{\pi}{4}$. The latter results in a solitary wave that is not aligned with the grid, so reflections occur at the boundaries of the domain immediately. For the simulation shown below, we have taken $u_0(x)=u_{\mathrm{kink}}(x,t_0)$, $\theta=\frac{2}{3}$, $k=20h$, $t_0=0$ and $t_f=20$. Moreover, we had to pick $\theta=\frac{2}{3}$ because for any $\theta\le\frac{1}{2}$ oscillations arose at the boundary, which are likely due to the scheme and not the equation, thus picking a value of $\theta$ a good bit into the "exponentially damped" spectrum of the time stepping schemes assures these oscillations are not created. 
+
+\image html christov.2d-angled_kink.gif "Animation of a moving 2D kink, at 45 degrees to the axis of the grid, showing boundary effects." width=5cm
+
+Another interesting solution to the sine-Gordon equation (which cannot be obtained analytically) can be produced by using two 1D breathers to construct the following separable 2D initial condition:
+\f[
+u_0(x) = u_{\mathrm{pseudobreather}}(x,t_0) = 16\arctan \left(\frac{m}{\sqrt{1-m^2}} \frac{\sin\left(\sqrt{1-m^2}t_0\right)}{\cosh(mx_1)} \right)\arctan \left(\frac{m}{\sqrt{1-m^2}} \frac{\sin\left(\sqrt{1-m^2}t_0\right)}{\cosh(mx_2)} \right),
+\f]
+where $x=(x_1,x_2)\in{R}^2$, $m=0.5(<1)$ as in the 1D case we discussed above. For the simulation shown below, we have chosen $\theta=\frac{1}{2}$, $k=10h$, $t_0=-5.4414$ and $t_f=2.7207$. The solution is pretty interesting --- it acts like a breather (as far as the pictures are concerned); however, it appears to break up and reassemble, rather than just oscillate. Naturally, the are things in the animation that can not be explained in a few words due to the fact that this initial condition is not an analytical solution to the sine-Gordon equation.
+
+\image html christov.2d-pseudobreather.gif "Animation of a 2D pseudobreather." width=5cm
+
+<h3>Future work</h3>
+Time-adaptive grids would likely be useful when solving a problem like this one (especially in the case of a traveling wave). However, due to the complexity of their implementation they are better left for another day.
diff --git a/deal.II/examples/step-25/step-25 b/deal.II/examples/step-25/step-25
new file mode 100755 (executable)
index 0000000..8b8ee8f
Binary files /dev/null and b/deal.II/examples/step-25/step-25 differ
diff --git a/deal.II/examples/step-25/step-25.cc b/deal.II/examples/step-25/step-25.cc
new file mode 100644 (file)
index 0000000..d3f1b9c
--- /dev/null
@@ -0,0 +1,819 @@
+/*    $Id: project.cc descends from heat-equation.cc, which       */
+/*         descends from step-4.cc (2006/03/01).                  */
+/*    Author: Ivan Christov, Texas A&M University, 2006           */
+
+/*    $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $        */
+/*    Version: $Name:  $                                          */
+/*                                                                */
+/*    Copyright (C) 2006 by the deal.II authors                   */
+/*                                                                */
+/*    This file is subject to QPL and may not be  distributed     */
+/*    without copyright and license information. Please refer     */
+/*    to the file deal.II/doc/license.html for the text and       */
+/*    further information on this license.                        */
+
+
+ // @sect3{Include files and global variables}
+
+// For an explanation of the include files, the reader should refer to
+// the example programs step-1 through step-4. They are in the
+// standard order, which is ``base'' -- ``lac'' -- ``grid'' --
+// ``dofs'' -- ``fe'' -- ``numerics'' (since each of these categories
+// roughly builds upon previous ones), then a few C++ headers for
+// file, input/output and string streams.
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_gmres.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.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/data_out_stack.h>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+// The following global variable is used to determine whether the
+// problem being solved is one for which an exact solution is known,
+// e.g. we are using the exact solution as the initial condition. It
+// is set to zero by default, and modified by ``InitialValues::value''
+// (see below). Things such as the computation of the error between
+// the numerical and exact solutions depend on the value of this
+// variable.
+bool exact_solution_known = false;
+
+// @sect3{The ``SineGordonProblem'' class template}
+
+// The entire algorithm for solving the problem is encapsulated in
+// this class. Also, note that the class is declared with a template
+// parameter, which is the spatial dimension, so that we can solve the
+// sine-Gordon equation in one, two or three spatial dimension. For
+// more on the dimension-independent class-encapsulation of the
+// problem, the reader should consult step-3 and step-4.
+template <int dim>
+class SineGordonProblem 
+{
+  public:
+    SineGordonProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void compute_nl_term (const Vector<double> &old_data, 
+                         const Vector<double> &new_data,
+                         Vector<double>       &nl_term) const;
+    void compute_nl_matrix (const Vector<double> &old_data, 
+                           const Vector<double> &new_data,
+                           SparseMatrix<double> &nl_matrix) const;
+    void solve ();
+    void compute_error  (const unsigned int timestep_number);
+    void output_results (const unsigned int timestep_number);
+
+    Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    SparseMatrix<double> mass_matrix;
+    SparseMatrix<double> laplace_matrix;
+    
+    double time, final_time, time_step;
+    double theta;
+
+    Vector<double>       solution, d_solution, old_solution;
+    Vector<double>       massmatxvel;
+    Vector<double>       system_rhs;
+    Vector<double>       fem_errors;
+
+    DataOutStack<dim>    data_out_stack;
+
+    static const unsigned int output_timestep_skip = 1;
+    static const int n_global_refinements = 6;
+};
+
+// @sect3{Exact solitary wave solutions of the sine-Gordon equation}
+
+// A kink-like solitary wave solution to the (``dim''+1) dimensional
+// sine-Gordon equation, which we can test our code against, is given
+// by Leibbrandt in \e Phys. \e Rev. \e Lett. \b 41(7), and is
+// implemented in the ``ExactSolution'' class.  However, it should be
+// noted that a closed-form solution can only be obtained for the
+// infinite-line initial-value problem (not the Neumann
+// initial-boundary-value problem under consideration here). However,
+// given that we impose \e zero Neumann boundary conditions, we expect
+// that the solution to our initial-boundary-value problem would be
+// close (in fact, equal) to the solution infinite-line initial-value
+// problem, if reflections of waves off the boundaries of our domain
+// do \e not occur.
+// 
+// The constants $\vartheta$ (``th'') and $\lambda$ (``ld'') in the 2D
+// solution and $\vartheta$ (``th''), $\phi$ (``phi'') and $\tau$
+// (``tau'') in the 3D solution are called the B&auml;cklund
+// transformation parameters. They control such things as the
+// orientation and steepness of the kink. For the purposes of testing
+// the code against the exact solution, one should choose the
+// parameters so that the kink is aligned with the grid, e.g. $\vartheta
+// = \phi = \pi$.
+//
+// In 1D, more interesting analytical solutions are known. Many of
+// them are listed on
+// http://mathworld.wolfram.com/Sine-GordonEquation.html . We have
+// implemented the one kink, two kink, kink-antikink and stationary
+// breather solitary-wave solutions.
+template <int dim>
+class ExactSolution : public Function<dim>
+{
+  public:
+    ExactSolution (const unsigned int n_components = 1,
+                  const double time = 0.) : Function<dim>(n_components, time) {};
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component = 0) const;
+};
+
+template <int dim>
+double ExactSolution<dim>::value (const Point<dim> &p,
+                                 const unsigned int /*component*/) const
+{
+  double t = this->get_time ();
+
+  if (dim==1)
+    {
+      double m = 0.5;
+//      double beta = std::sqrt(m*m-1.)/m;
+      double c1 = 0.;
+      double c2 = 0.;
+//       double s1 = 1.;
+//       double s2 = -1.;
+
+      /* one kink (m>1) */
+      /* return 4.*std::atan(std::exp(s1*(p[0]+s2*beta*t)/std::sqrt(1.-beta*beta))); */
+
+      /* two kinks (m>1) */
+      /* return 4.*std::atan(beta*std::sinh(beta*m*p[0])/std::cosh(beta*m*t)); */
+
+      /* kink-antikink (m>1) */
+      /* return -4.*std::atan(m/std::sqrt(m*m-1)*std::sinh(std::sqrt(m*m-1.)*t+c2)/
+        std::cosh(m*p[0]+c1)); */
+
+      /* stationary breather (m<1), period = 2.*pi*sqrt(1.-m*m) 
+        for m=0.5, -5.4414 <= t <= 2.7207 is a good time interval */
+      return -4.*std::atan(m/std::sqrt(1.-m*m)*std::sin(std::sqrt(1.-m*m)*t+c2)
+                          /std::cosh(m*p[0]+c1));
+    }
+  else if (dim==2)
+    {
+      double th  = deal_II_numbers::PI/4.;
+      double ld  = 1.;
+      double a0  = 1.;
+      double s   = 1.;
+      double arg = 0.;
+      arg = p[0]*std::cos(th) + std::sin(th)*(p[1]*std::cosh(ld)+t*std::sinh(ld));
+      return 4.*std::atan(a0*std::exp(s*arg));
+    }
+  else if (dim==3)
+    {
+      double th  = deal_II_numbers::PI;
+      double phi = deal_II_numbers::PI;
+      double tau = 1.;
+      double c0  = 1.;
+      double s   = 1.;
+      double arg = 0.;
+      arg = (p[0]*std::cos(th) + p[1]*std::sin(th)*std::cos(phi)
+            + std::sin(th)*std::sin(phi)*(p[2]*std::cosh(tau)+t*std::sinh(tau)));
+      return 4.*std::atan(c0*std::exp(s*arg));
+    }
+}
+
+// @sect3{Boundary values and initial values}
+
+// For our problem, we do not enforce Dirichlet boundary conditions
+// and the Neumann boundary conditions are enforced directly through
+// the variational formulation. However, since our problem is time
+// dependent, we must specify the value of the independent variable
+// $u$ at the initial time $t_0$. We do so via the ``InitialValues''
+// class below.
+template <int dim>
+class InitialValues : public Function<dim>
+{
+  public:
+    InitialValues (const unsigned int n_components = 1, 
+                  const double time = 0.) : Function<dim>(n_components, time) {};
+  
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component = 0) const;
+};
+
+template <int dim>
+double InitialValues<dim>::value (const Point<dim> &p,
+                                 const unsigned int /*component*/) const 
+{   
+  // We could also use a localized wave form for our initial
+  // condition, and see how it evolves when governed by the
+  // sine-Gordon equation. An example of such an initial condition is
+  // the following:
+  /*
+  exact_solution_known = false;
+  if ((p[0]>=-M_PI) && (p[0]<=M_PI) && (p[1]>=-M_PI) && (p[1]<=M_PI)) {
+    return std::cos(p[0]/2.)*std::cos(p[1]/2.);
+  } else {
+    return 0.;
+  }
+  */
+
+  // In 2D, another possibility for a localized-wave initial condition
+  // is a separable solution composed of two 1D breathers:
+  exact_solution_known = false;
+  double m = 0.5;
+  double t = this->get_time();
+  double argx = m/std::sqrt(1-m*m)*std::sin(std::sqrt(1-m*m)*t)/std::cosh(m*p[0]);
+  double argy = m/std::sqrt(1-m*m)*std::sin(std::sqrt(1-m*m)*t)/std::cosh(m*p[1]);
+  return 16.*std::atan(argx)*std::atan(argy);
+
+  // For the purposes of validating the program, we can use an exact
+  // solution of the sine-Gordon equation, at $t=t_0$, as the initial
+  // condition for our problem. Though, perhaps, this is not the most
+  // efficient way to implement the exact solution as the initial
+  // conditons, it is instuctive.
+  /*
+  exact_solution_known = false;
+  ExactSolution<dim> exact_solution (1, this->get_time());
+  return exact_solution.value (p);
+  */
+}
+
+// @sect3{Implementation of the ``SineGordonProblem'' class}
+
+// \b TO \b DO: present the big picture here?
+
+// @sect4{SineGordonProblem::SineGordonProblem}
+
+// This is the constructor of the ``SineGordonProblem'' class. It
+// specifies the desired polynomial degree of the finite elements,
+// associates a ``DoFHandler'' to the ``triangulation'' object (just
+// as in the example programs step-3 and step-4), initializes the
+// current or initial time, the final time, the time step size, and
+// the value of $\theta$ for the time stepping scheme.
+//
+// Note that if we were to chose the explicit Euler time stepping
+// scheme ($\theta = 0$), then we must pick a time step $k \le h$,
+// otherwise the scheme is not stable and oscillations might arise in
+// the solution. The Crank-Nicolson scheme ($\theta = \frac{1}{2}$)
+// and the implicit Euler scheme ($\theta=1$) do not suffer from this
+// deficiency, since they are unconditionally stable. However, even
+// then the time step should be chosen to be on the order of $h$ in
+// order to obtain a good solution.
+template <int dim>
+SineGordonProblem<dim>::SineGordonProblem () :
+                fe (1),
+               dof_handler (triangulation),
+               time (-5.4414/*0.*/),
+               final_time (2.7207/*20.*/),
+               time_step (10*1./std::pow(2.,n_global_refinements)),
+               theta (0.5)
+{}
+
+// @sect4{SineGordonProblem::make_grid_and_dofs}
+
+// This function creates a rectangular grid in ``dim'' dimensions and
+// refines it several times. Also, all matrix and vector members of
+// the ``SineGordonProblem'' class are initialized to their
+// approrpiate sizes once the degrees of freedom have been
+// assembled. Unlike its analogue in step-3 (and step-4) this function
+// uses ``MatrixCreator'' class to generate a mass matrix $M$ and a
+// Laplace matrix $A$ and store them in the appropriate variables
+// for the remainder of the program's life.
+template <int dim>
+void SineGordonProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, -10, 10);
+  triangulation.refine_global (n_global_refinements);
+  
+  std::cout << "   Number of active cells: "
+           << triangulation.n_active_cells()
+           << std::endl
+           << "   Total number of cells: "
+           << triangulation.n_cells()
+           << std::endl;
+
+  dof_handler.distribute_dofs (fe);
+
+  std::cout << "   Number of degrees of freedom: "
+           << dof_handler.n_dofs()
+           << std::endl;
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(), 
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress ();
+
+  system_matrix.reinit  (sparsity_pattern);
+  mass_matrix.reinit    (sparsity_pattern);
+  laplace_matrix.reinit (sparsity_pattern);
+
+  MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3), 
+                                    mass_matrix);
+  MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3), 
+                                       laplace_matrix);
+
+  solution.reinit       (dof_handler.n_dofs());
+  d_solution.reinit     (dof_handler.n_dofs());
+  old_solution.reinit   (dof_handler.n_dofs());
+  massmatxvel.reinit    (dof_handler.n_dofs());
+  system_rhs.reinit     (dof_handler.n_dofs());
+
+  // We will use the ``fem_errors'' vector, which is of size equal to
+  // the number of time steps, to store the errors in the finite
+  // element solution after each time step. Note that we must make the
+  // first element of the vector equal to zero, since there is no
+  // error in the solution after zeroth time step because the solution
+  // is just the initial condition.
+  const unsigned int n_time_steps
+    = static_cast<unsigned int>(std::ceil(std::fabs(final_time-time)/time_step));
+  fem_errors.reinit (n_time_steps);
+  fem_errors(0) = 0.;
+}
+
+// @sect4{SineGordonProblem::assemble_system}
+
+// This functions assembles the system matrix and right-hand side
+// vector for each iteration of Newton's method. The reader should
+// refer to the last section of the Introduction for the explicit
+// formulas for the system matrix and right-hand side.
+template <int dim>
+void SineGordonProblem<dim>::assemble_system () 
+{  
+  // First we assemble the Jacobian matrix $F'_h(U^n_l)$, where
+  // $U^n_l$ is stored in the vector ``solution'' for convenience.
+  system_matrix = 0;
+  system_matrix.copy_from (mass_matrix);
+  system_matrix.add (std::pow(time_step*theta,2), laplace_matrix);
+  SparseMatrix<double> tmp_matrix (sparsity_pattern);
+  compute_nl_matrix (old_solution, solution, tmp_matrix);
+  system_matrix.add (-std::pow(time_step*theta,2), tmp_matrix);
+
+  // Then, we compute the right-hand side vector $-F_h(U^n_l)$.
+  system_rhs = 0;
+
+  tmp_matrix = 0;
+  tmp_matrix.copy_from (mass_matrix);
+  tmp_matrix.add (std::pow(time_step*theta,2), laplace_matrix);
+  Vector<double> tmp_vector (solution.size());
+  tmp_matrix.vmult (tmp_vector, solution);
+  system_rhs += tmp_vector;
+
+  tmp_matrix = 0;
+  tmp_matrix.copy_from (mass_matrix);
+  tmp_matrix.add (-std::pow(time_step,2)*theta*(1-theta), laplace_matrix);
+  tmp_vector = 0;
+  tmp_matrix.vmult (tmp_vector, old_solution);
+  system_rhs -= tmp_vector;
+
+  system_rhs.add (-time_step, massmatxvel);
+
+  tmp_vector = 0;
+  compute_nl_term (old_solution, solution, tmp_vector);
+  system_rhs.add (std::pow(time_step,2)*theta, tmp_vector);
+
+  system_rhs *= -1;
+}
+
+// @sect4{SineGordonProblem::compute_nl_term}
+
+// This function computes the vector $S(\cdot,\cdot)$ corresponding to the
+// nonlinear term in the auxilliary (second) equation of the split
+// formulation. This function not only simplifies the repeated
+// computation of this term, but it is also a fundamental part of
+// nonlinear iterative solver that we use when the time stepping is
+// implicit (i.e. $\theta\ne 0$). Moreover, we must allow the function
+// to receive as input an "old" and a "new" solution, which may not be
+// the actual solutions of the problem stored in ``old_solution'' and
+// ``solution.'' For the purposes of this function, let us call the
+// first two arguments $w_{\mathrm{old}}$ and $w_{\mathrm{new}}$,
+// respectively.
+//
+// It is perhaps worth investigating what order quadrature formula is
+// best suited for this type of integration, since $\sin(\cdot)$ is an
+// oscillatory function.
+template <int dim>
+void SineGordonProblem<dim>::compute_nl_term (const Vector<double> &old_data,
+                                             const Vector<double> &new_data,
+                                             Vector<double>       &nl_term) const
+{
+  QGauss<dim>   quadrature_formula (3);
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values | update_JxW_values | update_q_points);
+  
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.n_quadrature_points;
+  
+  Vector<double> local_nl_term (dofs_per_cell);      
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell); 
+  std::vector<double> old_data_values (n_q_points);
+  std::vector<double> new_data_values (n_q_points);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    { 
+      // Once we re-initialize our ``FEValues'' instantiation to the
+      // current cell, we make use of the ``get_function_values''
+      // routine to get the obtain the values of the "old" data
+      // (presumably at $t=t_{n-1}$) and the "new" data (presumably at
+      // $t=t_n$) at the nodes of the chosen quadrature formula.
+      fe_values.reinit (cell);
+      fe_values.get_function_values (old_data, old_data_values);
+      fe_values.get_function_values (new_data, new_data_values);
+      
+      // Now, we can evaluate $\int_K \sin\left[\theta w_{\mathrm{new}} +
+      // (1-\theta) w_{\mathrm{old}}\right]\,\varphi_j\,\mathrm{d}x$ using
+      // the desired quadrature formula.
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)       
+         local_nl_term(i) += (std::sin(theta*new_data_values.at(q_point) +
+                                       (1-theta)*old_data_values.at(q_point)) *
+                              fe_values.shape_value (i, q_point) *
+                              fe_values.JxW (q_point));            
+      
+      // We conclude by adding up the contributions of the
+      // integrals over the cells to the global integral.
+      cell->get_dof_indices (local_dof_indices);
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       nl_term(local_dof_indices[i]) += local_nl_term(i);
+      
+      local_nl_term = 0;
+    }
+}
+
+// @sect4{SineGordonProblem::compute_nl_matrix}
+
+// This function computes the matrix $N(\cdot,\cdot)$ corresponding to
+// the nonlinear term in the Jacobian of $F(\cdot)$. It is also a
+// fundamental part of nonlinear iterative solver. Just as
+// ``compute_nl_term'', we must allow this function to receive
+// as input an "old" and a "new" solution, which we call the
+// $w_{\mathrm{old}}$ and $w_{\mathrm{new}}$, respectively.
+template <int dim>
+void SineGordonProblem<dim>::compute_nl_matrix (const Vector<double> &old_data, 
+                                               const Vector<double> &new_data,
+                                               SparseMatrix<double> &nl_matrix) const
+{
+  QGauss<dim>   quadrature_formula (3);
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values | update_JxW_values | update_q_points);
+  
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.n_quadrature_points;
+  
+  FullMatrix<double> local_nl_matrix (dofs_per_cell, dofs_per_cell);
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell); 
+  std::vector<double> old_data_values (n_q_points);
+  std::vector<double> new_data_values (n_q_points);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    { 
+      // Again, first we re-initialize our ``FEValues'' instantiation
+      // to the current cell.
+      fe_values.reinit (cell);
+      fe_values.get_function_values (old_data, old_data_values);
+      fe_values.get_function_values (new_data, new_data_values);
+      
+      // Then, we evaluate $\int_K \cos\left[\theta w_{\mathrm{new}} +
+      // (1-\theta) w_{\mathrm{old}}\right]\, \varphi_i\,
+      // \varphi_j\,\mathrm{d}x$ using the desired quadrature formula.    
+      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) 
+           local_nl_matrix(i,j) += (std::cos(theta*new_data_values.at(q_point) +
+                                             (1-theta)*old_data_values.at(q_point)) *
+                                    fe_values.shape_value (i, q_point) *
+                                    fe_values.shape_value (j, q_point) *
+                                    fe_values.JxW (q_point));
+      
+      // Finally, we add up the contributions of the integrals over
+      // the cells to the global integral.
+      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)
+         nl_matrix.add(local_dof_indices[i], local_dof_indices[j], 
+                       local_nl_matrix(i,j));
+
+      local_nl_matrix = 0;
+    }
+}
+
+// @sect4{SineGordonProblem::compute_error}
+
+// This function computes the norm of the difference between the
+// computed (i.e., finite element) solution after time step
+// ``timestep_number'' and the exact solution to see how well we are
+// doing. There are several choices for norms available to us in the
+// ``VectorTools'' class. We use the $L^2$ norm because it is a
+// natural choice for our problem, since the solutions to the
+// sine-Gordon equation have finite energy or, equivalently, are $L^2$
+// functions. Given our weak formulation of the sine-Gordon equation,
+// we are computing a solution $u\in H^1(\Omega)$, hence we could also
+// use the $H^1$ norm to compute the error of the spatial
+// discretization. For more information on the details behind this
+// computation, the reader should refer to step-7.
+template <int dim>
+void SineGordonProblem<dim>::compute_error (const unsigned int timestep_number)
+{
+  ExactSolution<dim> exact_solution (1, time);
+  
+  Vector<double> difference_per_cell (triangulation.n_active_cells());
+  VectorTools::integrate_difference (dof_handler,
+                                    solution,
+                                    exact_solution,
+                                    difference_per_cell,
+                                    QGauss<dim>(3),
+                                    VectorTools::L2_norm);
+  fem_errors(timestep_number) = difference_per_cell.l2_norm();
+  
+  std::cout << "   The L^2 error in the solution is " 
+           << fem_errors(timestep_number) << "."
+           << std::endl;
+}
+
+// @sect4{SineGordonProblem::solve}
+
+// This function uses the GMRES iterative solver on the linear system
+// of equations resulting from the finite element spatial
+// discretization of each iteration of Newton's method for the
+// (nonlinear) first equation in the split formulation we derived in
+// the Introduction. The solution to the system is, in fact, $\delta
+// U^n_l$ so it is stored in ``d_solution'' and used to update
+// ``solution'' in the ``run'' function. We cannot use the Conjugate
+// Gradient solver because the nonlinear term in the Jacobian matrix
+// results in a non-positive-definite matrix to invert. Moreover, we
+// would like the solver to quit when the \e relative error is
+// $10^{-12}$. This function is similar to its analogue in step-3 (and
+// step-4); the only difference is the choice of iterative solver and
+// the new stopping criterion.
+template <int dim>
+void SineGordonProblem<dim>::solve () 
+{
+  SolverControl solver_control (1000, 1e-12*system_rhs.l2_norm());
+  SolverGMRES<> gmres (solver_control);
+  d_solution = 0;
+  gmres.solve (system_matrix, d_solution, system_rhs, PreconditionIdentity());
+
+  std::cout << "   " << solver_control.last_step()
+           << " GMRES iterations needed to obtain convergence."
+           << std::endl;
+}
+
+// @sect4{SineGordonProblem::output_results}
+
+// This function outputs the results to a file. It is almost identical
+// to its counterpart in step-3 (and step-4). The only new thing is
+// that the function now takes a parameter --- the time step number
+// --- so that it can append it to the name of the file, which the
+// current solution is output to.
+template <int dim>
+void SineGordonProblem<dim>::output_results (const unsigned int timestep_number)
+{
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "u");
+  data_out.build_patches ();
+
+  std::ostringstream filename;
+  filename << "solution-" << dim << "d-";
+
+  // Pad the time step number in filename with zeros in the beginning
+  // so that the files are ordered correctly in the shell and we can
+  // generate a good animation using convert.
+  if (timestep_number<10) 
+    filename << "0000" << timestep_number;
+  else if (timestep_number<100)
+    filename << "000" << timestep_number;
+  else if (timestep_number<1000)
+    filename << "00" << timestep_number;
+  else if (timestep_number<10000)
+    filename << "0" << timestep_number;
+  else 
+    filename << timestep_number;
+  
+  // We output the solution at the desired times in ``vtk'' format, so
+  // that we can use VisIt to make plots and/or animations.
+  filename << ".vtk";
+  std::ofstream output (filename.str().c_str());
+  data_out.write_vtk (output);
+
+  // We also store the current solution in our instantiation of a
+  // ``DataOutStack'' object, so that we can make a space-time plot of
+  // the solution.
+  data_out_stack.new_parameter_value (time, time_step*output_timestep_skip);
+  data_out_stack.attach_dof_handler (dof_handler);
+  data_out_stack.add_data_vector (solution, "solution");
+  data_out_stack.build_patches (1);
+  data_out_stack.finish_parameter_value ();
+}
+
+// @sect4{SineGordonProblem::run}
+
+// This function has the top-level control over everything: it runs
+// the (outer) time-stepping loop, the (inner) nonlinear-solver loop,
+// outputs the solution after each time step and calls the
+// ``compute_error'' routine after each time step if an exact solution
+// is known.
+template <int dim>
+void SineGordonProblem<dim>::run () 
+{
+  data_out_stack.declare_data_vector ("solution",
+                                     DataOutStack<dim>::dof_vector);
+
+  std::cout << "Solving problem in " << dim << " space dimensions." 
+           << std::endl;
+  
+  make_grid_and_dofs ();
+
+  // To aknowledge the initial condition, we must use the function
+  // $u_0(x)$ to compute the zeroth time step solution $U^0$. Note
+  // that when we create the ``InitialValues'' ``Function'' object, we
+  // set its internal time variable to $t_0$, in case our initial
+  // condition is a function of space and time evaluated at $t=t_0$.
+  InitialValues<dim> initial_condition (1, time);
+
+  // Then, in 2D and 3D, we produce $U^0$ by projecting $u_0(x)$ onto
+  // the grid using ``VectorTools::project''. In 1D, however, we
+  // obtain the zeroth time step solution by interpolating $u_0(x)$ at
+  // the global degrees of freedom using
+  // ``VectorTools::interpolate''. We must make an exception for the
+  // 1D case because the projection algorithm computes integrals over
+  // the boundary of the domain, which do not make sense in 1D, so we
+  // cannot use it.
+  if (dim == 1) 
+    {
+      VectorTools::interpolate (dof_handler, initial_condition, solution);
+    }
+  else 
+    {
+      ConstraintMatrix constraints;
+      constraints.close();
+      VectorTools::project (dof_handler, constraints, QGauss<dim>(3),
+                           initial_condition, solution);
+    }
+
+  // For completeness, we output the zeroth time step to a file just
+  // like any other other time step.
+  output_results (0);
+
+  // Now we perform the time stepping: at every time step we solve the
+  // matrix equation(s) corresponding to the finite element
+  // discretization of the problem, and then advance our solution
+  // according to the time stepping formulas we discussed in the
+  // Introduction.
+  unsigned int timestep_number = 1;
+  for (time+=time_step; time<=final_time; time+=time_step, ++timestep_number)
+    {
+      old_solution = solution;
+
+      std::cout << std::endl
+               << " Time step #" << timestep_number << "; "
+               << "advancing to t = " << time << "." 
+               << std::endl;
+
+      // First we must solve the nonlinear equation in the split
+      // formulation via Newton's method --- i.e. solve for $\delta
+      // U^n_l$ then compute $U^n_{l+1}$ and so on. The stopping
+      // criterion is that $\|F_h(U^n_l)\|_2 \le 10^{-6}
+      // \|F_h(U^n_0)\|_2$. When the loop below is done, we have (an
+      // approximation of) $U^n$.
+      double initial_rhs_norm = 0.;
+      unsigned int nliter = 1;
+      do 
+       {
+         assemble_system ();
+         if (nliter == 1) initial_rhs_norm = system_rhs.l2_norm();
+         std::cout << "   [NLITER]"; 
+         solve ();
+         solution += d_solution;
+         d_solution = 0;
+         nliter++;
+       } 
+      while (system_rhs.l2_norm() > 1e-6 * initial_rhs_norm);
+  
+      // In the case of the explicit Euler time stepping scheme, we
+      // must pick the time step to be quite small in order for the
+      // scheme to be stable. Therefore, there are a lot of time steps
+      // during which "nothing interesting happens" in the
+      // solution. To improve overall efficiency --- in particular,
+      // speed up the program and save disk space --- we only output
+      // the solution after ``output_timestep_skip'' time steps have
+      // been taken.
+      if (timestep_number % output_timestep_skip == 0)
+         output_results (timestep_number);
+      
+      // Upon obtaining the solution to the problem at $t=t_n$, we
+      // must update the auxilliary velocity variable $V^n$. However,
+      // we do not compute and store $V^n$ since it is not a quantity
+      // we use directly in the problem. Hence, for simplicity, we
+      // update $MV^n$ directly using the second equation in the last
+      // subsection of the Introduction.
+      Vector<double> tmp_vector (solution.size());
+      laplace_matrix.vmult (tmp_vector, solution);
+      massmatxvel.add (-time_step*theta, tmp_vector);
+
+      tmp_vector = 0;
+      laplace_matrix.vmult (tmp_vector, old_solution);
+      massmatxvel.add (-time_step*(1-theta), tmp_vector);
+      
+      tmp_vector = 0;
+      compute_nl_term (old_solution, solution, tmp_vector);
+      massmatxvel.add (-time_step, tmp_vector);
+
+      // Before concluding the $n^{\mathrm{th}}$ time step, we compute
+      // the error in the finite element solution at $t=t_n$ if the
+      // exact solution to the problem being solved is known.
+      if (exact_solution_known)
+       compute_error (timestep_number);
+    }
+
+  // After the time stepping is complete, we report the maximum (over
+  // all time steps) of the errors in the finite element solution if
+  // the exact solution of the problem being solved is know.
+  if (exact_solution_known)
+    std::cout << " The maximum L^2 error in the solution was "
+             << fem_errors.linfty_norm() << "."
+             << std::endl << std::endl;
+
+  // Finally, we output the sequence of solutions stored
+  // ``data_out_stack'' to a file of the appropriate format.
+  std::ostringstream filename;
+  filename << "solution-" << dim << "d-" << "stacked" << ".vtk";
+  std::ofstream output (filename.str().c_str());
+  data_out_stack.write_vtk (output);
+}
+
+// @sect3{The ``main'' function}
+
+// This is the main function of the program. It creates an object of
+// top-level class and calls its principal function. Also, we supress
+// some of the library output by setting ``deallog.depth_console'' to
+// zero. Furthermore, if exceptions are thrown during the execution of
+// the run method of the ``SineGordonProblem'' class, we catch and
+// report them here. For more information about exceptions the reader
+// should consult step-6.
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      SineGordonProblem<2> sg_problem;
+      sg_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.