]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Removed unfinished step-41
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 11 Oct 2011 11:53:12 +0000 (11:53 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 11 Oct 2011 11:53:12 +0000 (11:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@24591 0785d39b-7218-0410-832d-ea1e28bc413d

62 files changed:
deal.II/examples/step-1/Makefile
deal.II/examples/step-34/parameters.prm
deal.II/examples/step-41/CamClay.prm [deleted file]
deal.II/examples/step-41/Makefile [deleted file]
deal.II/examples/step-41/doc/builds-on [deleted file]
deal.II/examples/step-41/doc/intro.dox [deleted file]
deal.II/examples/step-41/doc/kind [deleted file]
deal.II/examples/step-41/doc/results.dox [deleted file]
deal.II/examples/step-41/doc/tooltip [deleted file]
deal.II/examples/step-41/include/assemble_flags.h [deleted file]
deal.II/examples/step-41/include/base.h [deleted file]
deal.II/examples/step-41/include/base.templates.h [deleted file]
deal.II/examples/step-41/include/camclay.h [deleted file]
deal.II/examples/step-41/include/camclay.templates.h [deleted file]
deal.II/examples/step-41/include/domain.h [deleted file]
deal.II/examples/step-41/include/domain.templates.h [deleted file]
deal.II/examples/step-41/include/error_handler.h [deleted file]
deal.II/examples/step-41/include/error_handler.templates.h [deleted file]
deal.II/examples/step-41/include/linear_elastic.h [deleted file]
deal.II/examples/step-41/include/linear_elastic.templates.h [deleted file]
deal.II/examples/step-41/include/local_assemble.h [deleted file]
deal.II/examples/step-41/include/local_assemble.templates.h [deleted file]
deal.II/examples/step-41/include/local_assemble_base.h [deleted file]
deal.II/examples/step-41/include/local_assemble_base.templates.h [deleted file]
deal.II/examples/step-41/include/local_assemble_elastic_matrix.h [deleted file]
deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h [deleted file]
deal.II/examples/step-41/include/local_assemble_elastic_rhs.h [deleted file]
deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h [deleted file]
deal.II/examples/step-41/include/local_assemble_hardening_project.h [deleted file]
deal.II/examples/step-41/include/local_assemble_mass.h [deleted file]
deal.II/examples/step-41/include/local_assemble_mass.templates.h [deleted file]
deal.II/examples/step-41/include/local_assemble_plastic_project.h [deleted file]
deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h [deleted file]
deal.II/examples/step-41/include/local_assemble_scalar_project.h [deleted file]
deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h [deleted file]
deal.II/examples/step-41/include/my_tools.h [deleted file]
deal.II/examples/step-41/include/output_processor.h [deleted file]
deal.II/examples/step-41/include/output_processor.templates.h [deleted file]
deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h [deleted file]
deal.II/examples/step-41/include/parsed_symmetric_tensor_function.templates.h [deleted file]
deal.II/examples/step-41/include/utilities.h [deleted file]
deal.II/examples/step-41/include/vector_space.h [deleted file]
deal.II/examples/step-41/include/vector_space.templates.h [deleted file]
deal.II/examples/step-41/mesh/square.inp [deleted file]
deal.II/examples/step-41/source/base.cc [deleted file]
deal.II/examples/step-41/source/camclay.cc [deleted file]
deal.II/examples/step-41/source/domain.cc [deleted file]
deal.II/examples/step-41/source/error_handler.cc [deleted file]
deal.II/examples/step-41/source/linear_elastic.cc [deleted file]
deal.II/examples/step-41/source/local_assemble.cc [deleted file]
deal.II/examples/step-41/source/local_assemble_base.cc [deleted file]
deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc [deleted file]
deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc [deleted file]
deal.II/examples/step-41/source/local_assemble_mass.cc [deleted file]
deal.II/examples/step-41/source/local_assemble_plastic_project.cc [deleted file]
deal.II/examples/step-41/source/local_assemble_scalar_project.cc [deleted file]
deal.II/examples/step-41/source/main.cc [deleted file]
deal.II/examples/step-41/source/output_processor.cc [deleted file]
deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc [deleted file]
deal.II/examples/step-41/source/vector_space.cc [deleted file]
deal.II/examples/step-41/step-41.cc [deleted file]
deal.II/examples/step-6/step-6.cc

index 084f822e9a4696cb7883383e69d316827193e6a5..b4c21ed3257b7a625ac520999557efd7d26de80e 100644 (file)
@@ -33,8 +33,6 @@ D = ../../
 clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
 
 
-
-
 #
 #
 # Usually, you will not need to change anything beyond this point.
@@ -45,6 +43,41 @@ clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
 # settings
 include $D/common/Make.global_options
 
+CXX += -I$(VTK_DIR)/include/vtk-5.6
+LIBS += -L/usr/local/vtk/lib/vtk/ -I$(VTK_DIR)/include/vtk-5.6 \
+-lvtkCharts \
+-lvtkCommon \
+-lvtkDICOMParser \
+-lvtkFiltering \
+-lvtkGenericFiltering \
+-lvtkGeovis \
+-lvtkGraphics \
+-lvtkHybrid \
+-lvtkIO \
+-lvtkImaging \
+-lvtkInfovis \
+-lvtkNetCDF \
+-lvtkRendering \
+-lvtkViews \
+-lvtkVolumeRendering \
+-lvtkWidgets \
+-lvtkalglib \
+-lvtkexoIIc \
+-lvtkexpat \
+-lvtkfreetype \
+-lvtkftgl \
+-lvtkjpeg \
+-lvtklibxml2 \
+-lvtkmetaio \
+-lvtkpng \
+-lvtkproj4 \
+-lvtksqlite \
+-lvtksys \
+-lvtktiff \
+-lvtkverdict \
+-lvtkzlib 
+
+
 
 # Since the whole project consists of only one file, we need not
 # consider difficult dependencies. We only have to declare the
index 5e12caf422f9eb48ae23d0b519150e8248730cc3..03e374ec0a4f4239f2656aa0bd6572677f18eb71 100644 (file)
@@ -1,10 +1,10 @@
 # Listing of Parameters
 # ---------------------
 set Extend solution on the -2,2 box = true
-set External refinement             = 5
+set External refinement             = 7
 set Number of cycles                = 4
 set Run 2d simulation               = true
-set Run 3d simulation               = true
+set Run 3d simulation               = false
 
 
 subsection Exact solution 2d
diff --git a/deal.II/examples/step-41/CamClay.prm b/deal.II/examples/step-41/CamClay.prm
deleted file mode 100644 (file)
index 404fbd4..0000000
+++ /dev/null
@@ -1,306 +0,0 @@
-# Listing of Parameters
-# ---------------------
-set Compute error            = true
-set Error file format        = tex
-set Output error tables      = true
-
-# Comma separated list of names for the components. This will be used both for
-# error tables in text format and to output the solution to a file. Note that
-# in the case of a vector function the error name which is used to compute the
-# norm (supposing the type of the other components is 'Add') is the first
-# one.
-set Solution names           = u1,u2
-
-# Comma separated version of the same thing as above for the latex version of
-# the table.
-set Solution names for latex = u1,u2
-
-# Comma separated list of table names.
-set Table names              = error
-set Write error files        = false
-
-
-subsection Dirichlet Data
-  # Any constant used inside the function which is not a variable name.
-  set Function constants  = 
-
-  # Separate vector valued expressions by ';' as ',' is used internally by the
-  # function parser.
-  set Function expression = sin(2*pi*x)*cos(2*pi*y); cos(2*pi*x)*sin(2*pi*y)
-
-  # The name of the variables as they will be used in the function, separated
-  # by ','.
-  set Variable names      = x,y,t
-end
-
-
-subsection Domain Parameters
-  set Input mesh file            = square
-  set Input mesh format          = ucd
-  set Output mesh file           = square_out
-  set Path of domain mesh files  = mesh/
-
-  # If this is false, then the input mesh file below is ignored and a
-  # hyper-cube is created.
-  set Read domain mesh from file = false
-
-
-  subsection Grid Out Parameters
-    set Format = none
-
-
-    subsection DX
-      # Write all faces, not only boundary
-      set Write all faces = true
-
-      # Write the mesh connectivity as DX grid cells
-      set Write cells     = true
-
-      # If cells are written, additionally write their diameter as data for
-      # visualization
-      set Write diameter  = false
-
-      # Write faces of cells. These may be boundary faces or all faces between
-      # mesh cells, according to "Write all faces"
-      set Write faces     = false
-
-      # Write the volume of each cell as data
-      set Write measure   = false
-    end
-
-    subsection Eps
-      # Azimuth of the viw point, that is, the angle in the plane from the
-      # x-axis.
-      set Azimuth         = 30
-
-      # Number of points on boundary edges. Increase this beyond 2 to see
-      # curved boundaries.
-      set Boundary points = 2
-
-      # (2D only) Write cell numbers into the centers of cells
-      set Cell number     = false
-
-      # Draw lines with user flag set in different color
-      set Color by flag   = false
-
-      # Draw different colors according to grid level.
-      set Color by level  = false
-
-      # Elevation of the view point above the xy-plane.
-      set Elevation       = 30
-
-      # (2D only) if "Cell number" is true, writenumbers in the form
-      # level.number
-      set Level number    = false
-
-      # Width of the lines drawn in points
-      set Line width      = 0.5
-
-      # Size of the output in points
-      set Size            = 300
-
-      # Depending on this parameter, either thewidth or height of the eps is
-      # scaled to "Size"
-      set Size by         = width
-
-      # Write numbers for each vertex
-      set Vertex number   = false
-    end
-
-    subsection Gnuplot
-      set Boundary points = 2
-      set Cell number     = false
-    end
-
-    subsection Msh
-      set Write faces = false
-      set Write lines = false
-    end
-
-    subsection Ucd
-      set Write faces    = false
-      set Write lines    = false
-      set Write preamble = true
-    end
-
-    subsection XFig
-      set Boundary        = true
-      set Boundary points = 0
-      set Boundary style  = 0
-      set Boundary width  = 3
-      set Fill style      = 20
-      set Level color     = false
-      set Level depth     = true
-      set Line style      = 0
-      set Line width      = 1
-    end
-
-  end
-
-end
-
-
-subsection Elastic Moduli
-  # Any constant used inside the functions which is not a variable name.
-  set Function constants = 
-
-  # Separate different components expressions by ';' as ',' is used internally
-  # by the function parser.
-  set Row 1              = 0.5; 0.5; 0
-  set Row 2              = 0.5; 0
-  set Row 3              = 0.5
-
-  # The name of the variables as they will be used in the function, separated
-  # by ','.
-  set Variable names     = x,y,t
-end
-
-
-subsection Exact Solution
-  # Any constant used inside the function which is not a variable name.
-  set Function constants  = 
-
-  # Separate vector valued expressions by ';' as ',' is used internally by the
-  # function parser.
-  set Function expression = sin(2*pi*x)*cos(2*pi*y); cos(2*pi*x)*sin(2*pi*y)
-
-  # The name of the variables as they will be used in the function, separated
-  # by ','.
-  set Variable names      = x,y,t
-end
-
-
-subsection General Parameters
-  # Determines level of screen output
-  set Console Depth                      = 10
-
-  # Linear Solver Reduction Tolerance
-  set Linear Solver Reducation Tolerance = 1.0E-12
-
-  # Number of Convergence Cycles Past Initial Refinement
-  set Number of Convergence Cycles       = 3
-end
-
-
-subsection Loading Parameters
-  # End time for the simulation
-  set End Time for Simulation = 1.0
-
-  # Number of Loading Cycles to split time interval
-  set Number of Loading Steps = 3
-end
-
-
-subsection Neumann Data
-  # Any constant used inside the function which is not a variable name.
-  set Function constants  = 
-
-  # Separate vector valued expressions by ';' as ',' is used internally by the
-  # function parser.
-  set Function expression = 0; 0
-
-  # The name of the variables as they will be used in the function, separated
-  # by ','.
-  set Variable names      = x,y,t
-end
-
-
-subsection Table 0
-  # Evaluate convergence rates and add a column to the table for each computed
-  # norm.
-  set Add convergence rates          = true
-
-  # The extra columns to add to the table.
-  set Extra terms                    = cells,dofs
-
-  # The caption that will go under the table if we write the file in tex
-  # format. The default value for this object is the same name as the table
-  # name.
-  set Latex table caption            = error
-
-  # Each component is separated by a semicolon, and each norm by a comma.
-  # Implemented norms are Linfty, L2, W1infty, H1 and Add, which means that
-  # the norm is added to the previous component. Useful for vector valued
-  # functions.
-  set List of error norms to compute = Linfty, L2, H1 ; Linfty, L2, H1
-end
-
-
-subsection Vector Space Parameters
-  # Boundary indicator, followed by semicolomn and a list of components to
-  # which this boundary conditions apply. More boundary indicators can be
-  # separated by semicolumn. 1:0,1,4 ; 2,4:0,2
-  set Dirichlet boundary map = 0,1,2,3:0,1
-
-  # Ordering of the degrees of freedom: none, comp, cuth.
-  set Dof ordering           = cuth
-
-  # The finite element space to use. For vector finite elements use the
-  # notation FESystem[FE_Q(2)^2-FE_DGP(1)] (e.g. Navier-Stokes).
-  set Finite element space   = FESystem[FE_Q(1)^2]
-
-  # Degree of the mapping. If 0 is used, then a Cartesian mapping is assumed.
-  set Mapping degree         = 1
-
-  # Boundary indicators, followed by semicolomn and a list of components to
-  # which this boundary conditions apply. More boundary indicators can be
-  # separated by semicolumn. 1:0,1,4 ; 2,4:0,2
-  set Neumann boundary map   = 
-
-  # Boundary indicator, followed by semicolomn and a list of components to
-  # which this boundary conditions apply. More boundary indicators can be
-  # separated by semicolumn. 1:0,1,4 ; 2,4:0,2
-  set Other boundary map     = 3:0
-
-
-  subsection Grid Parameters
-    set Bottom fraction        = .3
-
-    # If this number is greater than zero, the mesh is distorted upon
-    # refinement in order to disrupt its structureness.
-    set Distortion coefficient = 0
-    set Global refinement      = 3
-
-    # A number of zero means no limit.
-    set Max number of cells    = 0
-
-    # fixed_number: the Top/Bottom threshold fraction of cells are flagged for
-    # refinement/coarsening. fixed_fraction: the cells whose error is
-    # Top/Bottom fraction of the total are refined/coarsened. optmized: try to
-    # reach optimal error distribution, assuming error is divided by 4 upon
-    # refining. global: refine all cells.
-    set Refinement strategy    = fixed_number
-    set Top fraction           = .3
-  end
-
-end
-
-
-subsection Yield Stress
-  # Any constant used inside the function which is not a variable name.
-  set Function constants  = 
-
-  # Separate vector valued expressions by ';' as ',' is used internally by the
-  # function parser.
-  set Function expression = 0
-
-  # The name of the variables as they will be used in the function, separated
-  # by ','.
-  set Variable names      = x,y,t
-end
-
-
-subsection f - Body Source
-  # Any constant used inside the function which is not a variable name.
-  set Function constants  = 
-
-  # Separate vector valued expressions by ';' as ',' is used internally by the
-  # function parser.
-  set Function expression = 8*pi*pi*sin(2*pi*x)*cos(2*pi*y); 8*pi*pi*cos(2*pi*x)*sin(2*pi*y)
-
-  # The name of the variables as they will be used in the function, separated
-  # by ','.
-  set Variable names      = x,y,t
-end
-
-
diff --git a/deal.II/examples/step-41/Makefile b/deal.II/examples/step-41/Makefile
deleted file mode 100644 (file)
index 6ea309f..0000000
+++ /dev/null
@@ -1,149 +0,0 @@
-# $Id$
-
-# The large projects Makefile looks much like the one for small
-# projects. Basically, only the following seven parameters need to be
-# set by you:
-
-application-name = camclay
-
-# The first denotes the dimension for which the program is to be
-# compiled:
-deal_II_dimension = 2
-
-# The second tells us the name of the executable. It is prefixed by
-# `lib/' to designate its destination directory. Note that the program
-# name depends on the dimension, so you can keep copies for the
-# different dimensions around:
-target   = lib/$(application-name)-$(deal_II_dimension)d
-
-# The `debug-mode' variable works as in the small projects Makefile:
-debug-mode = on
-
-# And so does the following variable. You will have to set it to
-# something more reasonable, of course.
-clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
-
-# Finally, here is a variable which tells the `run' rule which
-# parameters to pass to the executable. Usually, this will be the name
-# of an input file.
-run-parameters  = parameter-file.prm
-
-# Now, this is the last variable you need to set, namely the path to
-# the deal.II toplevel directory:
-D = ../../
-
-#
-#
-# Usually, you will not need to change anything beyond this point.
-#
-#
-# This tells `make' where to find the global settings and rules:
-include $D/common/Make.global_options
-
-
-################################################################
-# This example program will only work if Trilinos is installed. If this
-# is not the case, then simply redefine the main targets to do nothing
-ifneq ($(USE_CONTRIB_TRILINOS),yes)
-default run clean:
-       @echo
-       @echo "==========================================================="
-       @echo "=   This program cannot be compiled without Trilinos. Make="
-       @echo "=   sure you have Trilinos installed and detected during  ="
-       @echo "=   configuration of deal.II                              ="
-       @echo "==========================================================="
-       @echo
-
-else
-#
-################################################################
-
-
-# First get a list of files belonging to the project. Include files
-# are expected in `include/', while implementation files are expected
-# in `source/'. Object files are placed into `lib/[123]d', using the
-# same base name as the `.cc' file.
-cc-files    = $(shell echo source/*.cc)
-o-files     = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.$(OBJEXT))
-go-files    = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.g.$(OBJEXT))
-h-files     = $(wildcard include/*.h)
-lib-h-files = $(shell echo $D/include/deal.II/*/*.h)
-
-# Since the whole project consists of only one file, we need not
-# consider difficult dependencies. We only have to declare the
-# libraries which we want to link to the object file. deal.II has two
-# libraries: one for the debug mode version of the
-# application and one for optimized mode.
-libs.g   := $(lib-deal2.g)
-libs.o   := $(lib-deal2.o)
-
-
-
-# Now use the information from above to define the set of libraries to
-# link with and the flags to be passed to the compiler:
-ifeq ($(debug-mode),on)
-  libraries = $(go-files) $(libs.g)
-  flags     = $(CXXFLAGS.g)
-else
-  libraries = $(o-files) $(libs.o)
-  flags     = $(CXXFLAGS.o)
-endif
-
-
-# Then augment the compiler flags by a specification of the dimension
-# for which the program shall be compiled:
-flags += -Ddeal_II_dimension=$(deal_II_dimension)
-
-
-# The following two rules define how to compile C++ files into object
-# files:
-lib/$(deal_II_dimension)d/%.g.$(OBJEXT) :
-       @echo =====$(application-name)=======$(deal_II_dimension)d====debug=====$(MT)== $(<F)
-       @$(CXX) $(flags) -c $< -o $@
-lib/$(deal_II_dimension)d/%.$(OBJEXT) :
-       @echo =====$(application-name)=======$(deal_II_dimension)d====optimized=$(MT)== $(<F)
-       @$(CXX) $(flags) -c $< -o $@
-
-
-
-# Next define how to link the executable
-$(target)$(EXEEXT) : $(libraries) Makefile
-       @echo =====$(application-name)=======$(deal_II_dimension)d==============$(MT)== Linking $(@F)
-       @$(CXX) -o $@ $(libraries) $(LIBS) $(LDFLAGS)
-
-
-
-# Rule how to run the program
-run: $(target)$(EXEEXT)
-       ./$(target)$(EXEEXT) $(run-parameters)
-
-
-# Rule how to clean up. This is split into several different rules to
-# allow for parallel execution of commands:
-clean: clean-lib clean-data
-       -rm -f *~ */*~ */*/*~ lib/Makefile.dep
-
-clean-lib:
-       -rm -f lib/?d/*.$(OBJEXT) lib/?d/*.g.$(OBJEXT) $(target)$(EXEEXT) lib/TAGS
-
-clean-data:
-       -rm -f $(clean-up-files)
-
-
-# Again tell `make' which rules are not meant to produce files:
-.PHONY: clean clean-data clean-lib run
-
-
-# Finally produce the list of dependencies. Note that this time, the
-# object files end up in directories of their own, so we have to
-# modify the output a bit. The file with the dependencies is put into
-# `lib/'.
-lib/Makefile.dep: $(cc-files) $(h-files) $(lib-h-files) Makefile
-       @echo =====$(application-name)=======$(deal_II_dimension)d================== Remaking $@
-       @$D/common/scripts/make_dependencies $(INCLUDE) -Blib $(cc-files) \
-        | $(PERL) -p -e 's!^lib/(.*):!lib/$(deal_II_dimension)d/$$1:!g;' \
-               > $@
-
-include lib/Makefile.dep
-
-endif  # USE_CONTRIB_TRILINOS
diff --git a/deal.II/examples/step-41/doc/builds-on b/deal.II/examples/step-41/doc/builds-on
deleted file mode 100644 (file)
index 8b13789..0000000
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/deal.II/examples/step-41/doc/intro.dox b/deal.II/examples/step-41/doc/intro.dox
deleted file mode 100644 (file)
index 446bf8b..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-Introduction about Cam Clay plasticity models.
diff --git a/deal.II/examples/step-41/doc/kind b/deal.II/examples/step-41/doc/kind
deleted file mode 100644 (file)
index 56e049c..0000000
+++ /dev/null
@@ -1 +0,0 @@
-solids
diff --git a/deal.II/examples/step-41/doc/results.dox b/deal.II/examples/step-41/doc/results.dox
deleted file mode 100644 (file)
index 6121fb6..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-<h1>Results</h1>
-
-Astonishing results here.
diff --git a/deal.II/examples/step-41/doc/tooltip b/deal.II/examples/step-41/doc/tooltip
deleted file mode 100644 (file)
index 03e9264..0000000
+++ /dev/null
@@ -1 +0,0 @@
-Solving Cam Clay plasticity model, getting the Nobel Prize.
diff --git a/deal.II/examples/step-41/include/assemble_flags.h b/deal.II/examples/step-41/include/assemble_flags.h
deleted file mode 100644 (file)
index 544f432..0000000
+++ /dev/null
@@ -1,122 +0,0 @@
-//---------------------------------------------------------------------------
-//    $Id: fe_update_flags.h,v 1.31 2005/10/24 04:33:03 guido Exp $
-//    Version: $Name:  $
-//
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//---------------------------------------------------------------------------
-#ifndef __deal2__assemble_flags_h
-#define __deal2__assemble_flags_h
-
-#include <deal.II/base/config.h>
-/**
- * The enum type given to the constructors of LocalAssembleBase objects, 
- * telling those objects which data to assemble on each mesh cell. 
- * When the GlobalAssembler calls the local one, it checks for each flag, 
- * and if it finds one, it assemble the corresponding object.
- *
- * By default, all flags are off, i.e. no procedure will be called.
- *
- * You can select more than one flag by concatenation
- * using the bitwise or operator|(AssembleFlags,AssembleFlags).
- */
-enum AssembleFlags
-  {
-    //! No update
-    assemble_default                      = 0,
-    //! Cell term. 
-    /**
-     * Assemble the cell term. This is
-     * usually needed, unless the matrix
-     * is only a flux matrix.
-     */
-    assemble_cell                       = 0x0001,
-    //! Assemble boundary term.
-    /**
-     * Calls the assemble_boundary_term for
-     * each boundary face.     
-     */
-    assemble_boundary                  = 0x0002,
-    //! Assemble face term.
-    /** Call the assemble_face_term for each face of each cell in the
-     * triangulation
-     */
-    assemble_face                      = 0x0004,
-    /** Assemble rhs cell term. Used in assemble_rhs method.*/
-    assemble_rhs_cell                 = 0x0008,
-    /** Assemble rhs boundary terms. */
-    assemble_rhs_boundary             = 0x0010
-  };
-
-
-
-
-
-/**
- * Global operator which returns an object in which all bits are set
- * which are either set in the first or the second argument. This
- * operator exists since if it did not then the result of the bit-or
- * <tt>operator |</tt> would be an integer which would in turn trigger
- * a compiler warning when we tried to assign it to an object of type
- * AssembleFlags.
- */
-inline
-AssembleFlags
-operator | (AssembleFlags f1, AssembleFlags f2)
-{
-  return static_cast<AssembleFlags> (
-                                    static_cast<unsigned int> (f1) |
-                                    static_cast<unsigned int> (f2));
-}
-
-
-
-
-/**
- * Global operator which sets the bits from the second argument also
- * in the first one.
- */
-inline
-AssembleFlags &
-operator |= (AssembleFlags &f1, AssembleFlags f2)
-{
-  f1 = f1 | f2;
-  return f1;
-}
-
-
-/**
- * Global operator which returns an object in which all bits are set
- * which are set in the first as well as the second argument. This
- * operator exists since if it did not then the result of the bit-and
- * <tt>operator &</tt> would be an integer which would in turn trigger
- * a compiler warning when we tried to assign it to an object of type
- * AssembleFlags.
- */
-inline
-AssembleFlags
-operator & (AssembleFlags f1, AssembleFlags f2)
-{
-  return static_cast<AssembleFlags> (
-                                    static_cast<unsigned int> (f1) &
-                                    static_cast<unsigned int> (f2));
-}
-
-
-/**
- * Global operator which clears all the bits in the first argument if
- * they are not also set in the second argument.
- */
-inline
-AssembleFlags &
-operator &= (AssembleFlags &f1, AssembleFlags f2)
-{
-  f1 = f1 & f2;
-  return f1;
-}
-#endif
diff --git a/deal.II/examples/step-41/include/base.h b/deal.II/examples/step-41/include/base.h
deleted file mode 100644 (file)
index 1a13228..0000000
+++ /dev/null
@@ -1,113 +0,0 @@
-#ifndef BASE_H
-#define BASE_H
-
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/base/parsed_function.h>
-#include <deal.II/base/timer.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <deal.II/grid/filtered_iterator.h>
-
-#include "vector_space.h"
-#include "domain.h"
-#include "error_handler.h"
-#include "linear_elastic.h"
-#include "camclay.h"
-
-
-
-using namespace dealii;
-using namespace dealii::Functions;
-
-template <int dim>
-class Base 
-{
- public:
-  Base ();
-  ~Base ();
-
-  //the routine called from int main()
-  void run ();
-  
- private:
-
-  //function with actual work, called from run()
-  void run_cc (const unsigned int &cc);
-
-  //function with actual work, called from run()
-  void run_step (const unsigned int &cc,
-                const unsigned int &step);
-  
-  //parse the parameters and get things rollin
-  void parameters ();
-
-  //elastic predictor
-  void elastic_predictor (double &time);
-
-  //plastic corrector
-  void plastic_corrector ();
-  
-  //writes the files
-  void write_files (const unsigned int &cc,
-                   const unsigned int &step,
-                   const Vector<double> &my_elas,
-                   const Vector<double> &my_plas,
-                   const Vector<double> &my_other);
-
-  void write_stress_strain (const unsigned int &cc,
-                           const unsigned int &step,
-                           const Vector<double> &my_elas);
-
-  void write_plot_values(const unsigned int &cc,
-                        const unsigned int &step,
-                        const Vector<double> &my_elas,
-                        const Vector<double> &my_hardening,
-                        const Vector<double> &my_plastic,
-                        const ParsedSymmetricTensorFunction<4,dim> &C);
-  
-  //A parameter handler object
-  ParameterHandler prm;
-
-  //ErrorHandler Object
-  ErrorHandler<dim, Vector<double> >  error_handler;
-  
-  //This holds the triangulation information
-  Domain<dim> domain;
-
-  //A vectorspace for the thermoelasticity
-  VectorSpace<dim> vspace;
-
-  //Elasticity Class - to change models, just change the class
-  LinearElastic<dim> elastic;
-  //HypoElastic<dim> elastic;
-  //HyperElastic<dim> elastic;
-  //LinearThermoElastic<dim> elastic;
-
-  //Plasticity Class - to change models, just change the class name
-  CamClay<dim> plastic;
-  //CamClayExplicit<dim> plastic;
-  //CamClayInvariant<dim> plastic;
-  //CamClayShearBand<dim> plastic;
-  //J2Flow<dim> plastic;
-  //J2FlowThermal<dim> plastic;
-  
-  //General Parameters
-  double lin_red_tol;
-  unsigned int num_steps;
-  double end_time;
-  unsigned int num_cc;
-  unsigned int num_threads;
-  unsigned int console_depth;
-
-  //various functions for the laplace equation
-  ParsedFunction<dim> exact_solution;
-
-
-};
-
-  
-#endif
diff --git a/deal.II/examples/step-41/include/base.templates.h b/deal.II/examples/step-41/include/base.templates.h
deleted file mode 100644 (file)
index 5f6f4c6..0000000
+++ /dev/null
@@ -1,468 +0,0 @@
-#include "../include/base.h"
-
-template <int dim>
-Base<dim>::Base ()
-  :
-  exact_solution(dim)
-{}
-
-template <int dim>
-Base<dim>::~Base () 
-{}
-
-template <int dim>
-void Base<dim>::run()
-{
-  //setup the depth of screen reporting
-  deallog.depth_console (10);
-
-  //get the parameters
-  parameters ();
-
-  //run as many convergence cycles as we need
-  for(unsigned int cc=0; cc<num_cc; ++cc) {
-    run_cc (cc);
-  }
-
-  //at the end of the problem, output the error table
-  error_handler.output_table(0);
-  
-}
-
-template <int dim>
-void Base<dim>::run_cc(const unsigned int &cc)
-{
-  deallog.push("CC-"+Utilities::int_to_string(cc));
-  deallog << "###############  Convergence Cycle: " 
-         << Utilities::int_to_string(cc) 
-         << " ###############" << std::endl;
-  
-  //for each convergence cycle, refine
-  //for now, we just leave things global
-  if(cc > 0) vspace.get_tria().refine_global(1);
-
-  vspace.redistribute_dofs();
-  
-  elastic.reinit(vspace);
-  
-  plastic.reinit(vspace);
-
-  plastic.initial_conditions(vspace);
-
-  deallog << "Building Matrix..." << std::endl;
-  elastic.build_matrix(prm, vspace);
-  plastic.build_matrix(vspace);
-
-  write_files (cc, 0,
-              elastic.sol_total,
-              plastic.solution,
-              plastic.sol_hard_iter);
-
-  write_plot_values(cc, 0, elastic.sol_total, plastic.sol_hard_iter,
-                   plastic.solution, elastic.C);
-
-  //just loop over the loading steps
-  for(unsigned int step=1; step<num_steps; ++step) {
-    
-    run_step(cc, step);
-    
-  }
-
-  //we are interested in the errors, after the whole run
-  exact_solution.set_time(end_time);
-  
-  error_handler.error_from_exact(vspace.get_dh(),
-                                elastic.sol_total,
-                                exact_solution,
-                                0, 0);
-  
-  deallog.pop();
-
-}
-
-template <int dim>
-void Base<dim>::run_step(const unsigned int &cc,
-                        const unsigned int &step)
-{
-  deallog << std::endl;
-  deallog.push("Step-" + Utilities::int_to_string(step,2));
-
-  double time = end_time * (double(step)/double(num_steps));
-
-  deallog << "Time: " << time <<std::endl;
-
-  elastic_predictor (time);
-
-  plastic_corrector ();
-
-  write_files (cc, step,
-              elastic.sol_total,
-              plastic.solution,
-              plastic.sol_hard_iter);
-
-  write_stress_strain (cc, step, elastic.sol_total);
-
-  write_plot_values(cc, step, elastic.sol_total, plastic.sol_hard_iter,
-                   plastic.solution, elastic.C);
-
-  deallog.pop();
-}
-
-template <int dim>
-void Base<dim>::elastic_predictor (double &time)
-{
-  deallog.push("ELASTIC");
-
-  elastic.reinit_step(time);
-
-  elastic.build_rhs(vspace);
-
-  elastic.solve(vspace, lin_red_tol);
-
-  deallog.pop();
-}
-
-template <int dim>
-void Base<dim>::plastic_corrector ()
-{
-  deallog.push("PLASTIC");
-  
-  plastic.reinit_step();
-  
-  plastic.update_internal_variables(vspace,
-                                   elastic.sol_total,
-                                    elastic.C);
-  
-  Timer timer;
-  timer.start();
-  plastic.project_strain(prm, vspace, elastic.A);
-  plastic.project_hardening(vspace);
-  timer.stop();
-  deallog << "Time for plastic projection: " << timer() << std::endl;
-  
-  deallog.pop();
-}
-
-
-template<int dim>
-void Base<dim>::write_files (const unsigned int &cc,
-                            const unsigned int &step,
-                            const Vector<double> &my_elas,
-                            const Vector<double> &my_plas,
-                            const Vector<double> &my_other)
-{
-  DataOut<dim> data_out;
-  data_out.attach_dof_handler (vspace.get_dh());
-
-  std::vector<std::string> elastic_names;
-  std::vector<std::string> plastic_names;
-  std::vector<std::string> other_names;
-  if(dim == 1) {
-    elastic_names.push_back("u1_elastic");
-    plastic_names.push_back("u1_plastic");
-  } 
-
-  if(dim == 2) {
-    elastic_names.push_back("u1_elastic");
-    elastic_names.push_back("u2_elastic");
-    plastic_names.push_back("u1_plastic");
-    plastic_names.push_back("u2_plastic");
-    other_names.push_back("hardening");
-    other_names.push_back("Newton_Iterations");
-  } 
-
-  if(dim == 3) {
-    elastic_names.push_back("u1_elastic");
-    elastic_names.push_back("u2_elastic");
-    elastic_names.push_back("u3_elastic");
-    plastic_names.push_back("u1_plastic");
-    plastic_names.push_back("u2_plastic");
-    plastic_names.push_back("u3_plastic");
-    other_names.push_back("hardening");
-    other_names.push_back("Newton_Iterations");
-    other_names.push_back("Nothing");
-  } 
-
-  data_out.add_data_vector(my_elas, elastic_names);
-  data_out.add_data_vector(my_plas, plastic_names);
-  data_out.add_data_vector(my_other, other_names);
-  data_out.add_data_vector(plastic.iterations, "Average_Newton_Iterations");
-
-  std::vector<unsigned int> partition_int (vspace.get_tria().n_active_cells());
-  GridTools::get_subdomain_association (vspace.get_tria(), partition_int);
-  const Vector<double> partitioning(partition_int.begin(),
-                                   partition_int.end());
-  data_out.add_data_vector (partitioning, "partitioning");
-
-  data_out.build_patches (vspace.get_fe().degree);
-
-  std::ostringstream filename;
-  filename << "out/solution-";
-  //first the convergence cycle
-  filename << std::setfill('0');
-  filename << std::setw(2) << cc;
-  filename << "-";
-  //then the step
-  filename << std::setfill('0');
-  filename << std::setw(5) << step;
-  filename << ".vtk";
-
-  std::ofstream output (filename.str().c_str());
-  data_out.write_vtk (output);
-
-}
-
-template<int dim>
-void Base<dim>::write_stress_strain (const unsigned int &cc,
-                                    const unsigned int &step,
-                                    const Vector<double> &my_elas)
-{
-  //reopen the just written file
-  std::ostringstream filename;
-  filename << "out/solution-";
-  //first the convergence cycle
-  filename << std::setfill('0');
-  filename << std::setw(2) << cc;
-  filename << "-";
-  //then the step
-  filename << std::setfill('0');
-  filename << std::setw(5) << step;
-  filename << ".vtk";
-  std::ofstream output (filename.str().c_str(),std::ios::app);
-
-  //WARNING - THE FOLLOWING IS NOT DIM INDEPENDENT!!!
-
-  //create an fe object to get the gradient values
-  QIterated<dim> qf(QTrapez<1>(),2);
-  FEValues<dim> fe_v(vspace.get_fe(), qf,
-                    update_gradients | update_quadrature_points);
-  std::vector< std::vector< Tensor<1,dim> > > total_grads(qf.size(),
-                                                         std::vector<Tensor<1,dim> >(dim));
-  std::vector< Point<dim> > points(qf.size());
-
-  output << "SCALARS uyy double 1" << std::endl;
-  output << "LOOKUP_TABLE default" << std::endl;
-
-  typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
-    endc = vspace.get_dh().end();
-  for (; cell!=endc; ++cell) {
-
-    fe_v.reinit(cell);
-    fe_v.get_function_gradients(my_elas, total_grads);
-
-    points = fe_v.get_quadrature_points();
-
-    //loop over the quadrature points---------------------------------
-    for (unsigned int qp = 0; qp<qf.size(); ++qp) {
-
-      output << total_grads[qp][1][1] << " ";
-
-    }
-  }
-
-
-  //be responsible - close your files
-  output.close();
-
-}
-
-template<int dim>
-void Base<dim>::write_plot_values(const unsigned int &cc,
-                                 const unsigned int &step,
-                                 const Vector<double> &my_elas,
-                                 const Vector<double> &my_hardening,
-                                 const Vector<double> &my_plas,
-                                 const ParsedSymmetricTensorFunction<4,dim> &C)
-{
-  double p=0;
-  double q=0;
-  double deve=0;
-  double k = 0;
-  //std::vector< std::vector< double> > s(dim, std::vector<double>(dim));
-
-  //reopen the just written file
-  std::ostringstream filename;
-  filename << "plot_data/data-";
-  //first the convergence cycle
-  filename << std::setfill('0');
-  filename << std::setw(2) << cc;
-  filename << ".txt";
-  std::ofstream output;
-
-  if(step == 0) {
-   output.open(filename.str().c_str(),std::ios::trunc);
-  } else {
-   output.open(filename.str().c_str(),std::ios::app);
-  }
-
-  //create an fe object to get the gradient values at the nodes
-  //QIterated<dim> qf(QTrapez<1>(),1);
-  QTrapez<dim> qf;
-  FEValues<dim> fe_v(vspace.get_fe(), qf,
-                    update_values | update_gradients | update_quadrature_points);
-  std::vector< std::vector< Tensor<1,dim> > > elastic_grads(qf.size(),
-                                                           std::vector<Tensor<1,dim> >(dim));
-  std::vector< std::vector< Tensor<1,dim> > > plastic_grads(qf.size(),
-                                                           std::vector<Tensor<1,dim> >(dim));
-  std::vector< Vector<double> > hard_values(qf.size(), Vector<double>(dim));
-  std::vector< Point<dim> > points(qf.size());
-
-  //std::vector< std::vector<double> > strain(dim, std::vector<double>(dim) );
-
-  //calculate p,q,deve
-  typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
-    endc = vspace.get_dh().end();
-  for (; cell!=endc; ++cell) {
-
-    fe_v.reinit(cell);
-
-    if (cell->index() == 20 ) {
-      fe_v.get_function_gradients(my_elas, elastic_grads);
-      fe_v.get_function_gradients(my_plas, plastic_grads);
-      fe_v.get_function_values(my_hardening, hard_values);
-      
-      points = fe_v.get_quadrature_points();
-
-      SymmetricTensor<2,dim> strain;
-      SymmetricTensor<2,dim> stress;
-      SymmetricTensor<2,3> s;
-
-      for(unsigned int a=0; a<dim; ++a) {
-       for(unsigned int b=0; b<dim; ++b) {
-         strain[a][b] = ( elastic_grads[0][a][b] - plastic_grads[0][a][b]);
-       }
-      }
-
-      for(unsigned int a=0; a<dim; ++a) {
-       for(unsigned int b=0; b<dim; ++b) {
-         for(unsigned int m=0; m<dim; ++m) {
-           for(unsigned int n=0; n<dim; ++n) {
-             stress[a][b] = C(points[0])[a][b][m][n] * strain[m][n];
-           }
-         }
-       }
-      }
-      //explicitly assuming the loading/symmetry condition
-
-      p = (-1.0/3.0) * (2*stress[0][0] + stress[1][1]);
-
-      s[0][0] = stress[0][0];
-      s[1][1] = stress[0][0];
-      s[2][2] = stress[1][1];
-
-      for(unsigned int kk=0; kk<3; ++kk) s[kk][kk] += p;
-
-      q = std::sqrt(3.0/2.0) * s.norm();
-
-      k = hard_values[0](0);
-
-    }
-  }
-
-  //write what we want
-  output << step << " "
-        << p << " "
-        << q << " "
-        << deve << " " 
-        << k << " " << std::endl;
-
-
-  //be responsible - close your files
-  output.close();
-
-}
-
-template<int dim>
-void Base<dim>::parameters ()
-{
-  deallog.push("PARAMETERS");
-
-  domain.declare_parameters(prm);
-
-  //add settings for the elastic vector space to the prm file
-  std::string space_name = "Vector Space Parameters";
-  vspace.declare_parameters(prm, space_name);
-
-  error_handler.declare_parameters(prm);
-
-  elastic.declare_parameters(prm);
-
-  plastic.declare_parameters(prm);
-  
-  prm.enter_subsection("General Parameters");
-  prm.declare_entry ("Linear Solver Reducation Tolerance",
-                    "1.0E-16",
-                    Patterns::Double(),
-                    "Linear Solver Reduction Tolerance");
-  prm.declare_entry ("Number of Convergence Cycles",
-                    "1",
-                    Patterns::Integer(),
-                    "Number of Convergence Cycles Past Initial Refinement");
-  prm.declare_entry ("Console Depth",
-                    "10",
-                    Patterns::Integer(),
-                    "Determines level of screen output");
-  prm.leave_subsection();
-
-  prm.enter_subsection("Loading Parameters");
-  prm.declare_entry ("End Time for Simulation",
-                    "1.0",
-                    Patterns::Double(),
-                    "End time for the simulation");
-  prm.declare_entry ("Number of Loading Steps",
-                    "10",
-                    Patterns::Integer(),
-                    "Number of Loading Cycles to split time interval");
-  prm.leave_subsection();
-
-  prm.enter_subsection("Exact Solution");
-  Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
-  prm.leave_subsection();
-  
-  // ==============================
-  
-  prm.read_input("CamClay.prm");
-  
-  // ==============================
-  //Initialize the domain    
-  domain.reinit(prm);
-  /*
-  domain.get_tria().clear();
-
-  Point<dim> lowerleft(0.0,0.0);
-  Point<dim> upperright(0.3,0.1);
-  std::vector<unsigned int> subs(2);
-  subs[0] = 30; subs[1]=10;
-  GridGenerator::subdivided_hyper_rectangle(domain.get_tria(),
-                                                subs,
-                                                lowerleft,
-                                                upperright, true);
-  */  
-  vspace.reinit(prm, domain.get_tria());
-  
-  error_handler.parse_parameters(prm);
-
-  elastic.parse_parameters(prm);
-
-  plastic.parse_parameters(prm);
-  
-  prm.enter_subsection("Exact Solution");
-  exact_solution.parse_parameters(prm);
-  prm.leave_subsection();
-  
-  prm.enter_subsection("General Parameters");
-  lin_red_tol = prm.get_double("Linear Solver Reducation Tolerance");
-  num_cc = prm.get_integer("Number of Convergence Cycles");
-  console_depth = prm.get_integer("Console Depth");
-  prm.leave_subsection();
-
-  prm.enter_subsection("Loading Parameters");
-  num_steps = prm.get_integer("Number of Loading Steps");
-  end_time = prm.get_double("End Time for Simulation");
-  prm.leave_subsection();
-
-  deallog.pop();
-
-}
-
-
diff --git a/deal.II/examples/step-41/include/camclay.h b/deal.II/examples/step-41/include/camclay.h
deleted file mode 100644 (file)
index e1ec253..0000000
+++ /dev/null
@@ -1,172 +0,0 @@
-#ifndef CAMCLAY
-#define CAMCLAY
-
-//deal.ii packages
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/base/parsed_function.h>
-#include <deal.II/base/symmetric_tensor.h>
-#include <deal.II/base/table.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/lac/sparse_direct.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/solver_control.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/numerics/matrices.h>
-
-//Trilinos packages
-#include <Sacado.hpp>
-
-//LPCM includes
-#include "vector_space.h"
-#include "local_assemble_plastic_project.h"
-#include "local_assemble_scalar_project.h"
-#include "my_tools.h"
-#include "domain.h"
-
-using namespace dealii;
-using namespace dealii::Functions;
-
-template <int dim>
-class CamClay
-{
- public:
-  
-  CamClay(); 
-  ~CamClay(); 
-       
-  void declare_parameters(ParameterHandler &prm);
-
-  void parse_parameters(ParameterHandler &prm);
-
-  void reinit(VectorSpace<dim> &vspace);
-
-  void reinit_step();
-
-  void project_strain(ParameterHandler &prm,
-                     VectorSpace<dim> &vspace,
-                     SparseMatrix<double> &A);
-
-  void project_hardening(VectorSpace<dim> &vspace);
-
-  void build_matrix(VectorSpace<dim> &vspace);
-
-  void initial_conditions(VectorSpace<dim> &vspace);
-
-  void update_internal_variables(VectorSpace<dim> &vspace,
-                                Vector<double> &elastic_solution,
-                                 ParsedSymmetricTensorFunction<4,dim> &C);
-
-  void compute_jacobian(FullMatrix<double> &jac,
-                        const Vector<double> &sol,
-                        const double p,
-                        const double q,
-                        const double M,
-                        const SymmetricTensor<4,dim> &C_qp,
-                        const SymmetricTensor<2,dim> &xi);
-
-  Vector<double> dF_dstress(const SymmetricTensor<2,dim> &stress,
-                            const double k);
-
-  //this is the plastic strain projected into the displacement  
-  Vector<double> solution;
-
-  //this is the average number of iterations per cell for the current step  
-  Vector<double> iterations;
-  
-  //internal variable vectors
-  Table<2, SymmetricTensor<2,dim> > plastic_strain;
-  Table<2, double > hardening;
-  Table<2, double> iter_table;
-
-  Vector<double> sol_hard_iter;
-               
- private:
-
-  //computes the yeild function
-  inline double yield_function(const double p, const double q,
-                              const double k);
-
-  //does the actual solving for the new internal variables
-  double solve(const int &index,
-               const unsigned int &qp,
-               const SymmetricTensor<2,dim> &trial_strain);
-  
-  //the potential for the hardening
-  double h(const double p,
-          const double k);
-  
-  //computes the hydrostatic stress invariant
-  inline double p(const SymmetricTensor<2,dim> &stress);
-
-  //computes the deviatoric stress invariant
-  inline double q(const SymmetricTensor<2,3> &xi);
-
-  //computes the deviatoric stress tensor
-  inline SymmetricTensor<2,3> xi(const SymmetricTensor<2,dim> &stress);
-
-  //computer_theresidua and the jacobian
-  void compute_res_jac(Vector<double> &res,
-                      FullMatrix<double> &jac,
-                      const Vector<double> &sol,
-                      const Vector<double> &prev_step,
-                      const SymmetricTensor<2,dim> &total_strain);
-  /*
-  //computes the residual
-  void compute_residual(Vector<double> &res,
-                        const Vector<double> &prev_step,
-                        const Vector<double> &sol,
-                        const SymmetricTensor<2,dim> &stress);
-
-
-
-  //creates the jacobian numerical
-  void numerical_jacobian(FullMatrix<double> &jac,
-                         const Vector<double> &sol,
-                         const Vector<double> &prev_sol,
-                         const SymmetricTensor<2,dim> &total_strain);*/
-
-  //translates from sym tensor notation to voigt notation
-  inline unsigned int sym2voigt(const unsigned int i,
-                               const unsigned int j);
-
-  inline unsigned int sac_num(const unsigned int i,
-                             const unsigned int j);
-
-  //translates voigt notation to sym tensor indices
-  inline std::vector<unsigned int> voigt2sym(const unsigned int i);
-
-  //delta tensor
-  inline double delta(const unsigned int i,
-                     const unsigned int j);
-
-  //voigt delta tensor
-  inline double delta(const unsigned int i);
-
-  //the yeild stress function - not needed for cam clay actually
-  ParsedFunction<dim> yield_stress;
-
-  //some constants
-  //double shear_mod;
-  //double bulk_mod;
-  double M;
-
-  SymmetricTensor<4,dim> C_qp;
-
-  SparseMatrix<double> MM;
-
-  SparsityPattern sp_MM;
-
-
-  
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/camclay.templates.h b/deal.II/examples/step-41/include/camclay.templates.h
deleted file mode 100644 (file)
index ca9c2ba..0000000
+++ /dev/null
@@ -1,670 +0,0 @@
-#include "../include/camclay.h"
-
-template <int dim>
-CamClay<dim>::CamClay()
- :
-yield_stress(1)
-{}
-
-template <int dim>
-CamClay<dim>::~CamClay()
-{
-  MM.clear();
-}
-
-template <int dim>
-void CamClay<dim>::declare_parameters(ParameterHandler &prm)
-{
-  prm.enter_subsection("Yield Stress");
-  Functions::ParsedFunction<dim>::declare_parameters(prm, 1);
-  prm.leave_subsection();
-
-}
-
-template <int dim>
-void CamClay<dim>::parse_parameters(ParameterHandler &prm)
-{
-  prm.enter_subsection("Yield Stress");
-  yield_stress.parse_parameters(prm);
-  prm.leave_subsection();
-}
-
-template <int dim>
-void CamClay<dim>::reinit(VectorSpace<dim> &vspace)
-{
-
-  //the total number of quadrature points
-  double total_qp = std::pow(double(2 * vspace.get_fe().degree + 1), dim);
-
-  //Resize the tables to hold the plastic strains and other internal variables.
-  plastic_strain.reinit(vspace.get_tria().n_active_cells(), int(total_qp));
-  hardening.reinit(vspace.get_tria().n_active_cells(), int(total_qp));
-  iter_table.reinit(vspace.get_tria().n_active_cells(), int(total_qp));
-
-  solution.reinit(vspace.n_dofs());
-  sol_hard_iter.reinit(vspace.n_dofs());
-
-  iterations.reinit(vspace.get_tria().n_active_cells());
-  iterations = 0; 
-
-  //setup for projection mass matrix
-  sp_MM.reinit(vspace.n_dofs(), 
-              vspace.n_dofs(), 
-              vspace.get_dh().max_couplings_between_dofs());
-
-  DoFTools::make_sparsity_pattern (static_cast<DoFHandler<dim> &> (vspace.get_dh()), sp_MM);
-  sp_MM.compress();
-
-  MM.reinit(sp_MM);
-}
-
-template <int dim>
-void CamClay<dim>::reinit_step()
-{
-
-  solution = 0;
-  iterations = 0;
-  sol_hard_iter = 0;
-
-}
-
-template <int dim>
-void CamClay<dim>::initial_conditions(VectorSpace<dim> &vspace)
-{
-
-  QGauss<dim>   qf_v((2*vspace.get_fe().degree) + 1);
-  
-  FEValues<dim> fe_v (vspace.get_fe(), qf_v, 
-                     update_values   | update_gradients |
-                     update_quadrature_points | update_JxW_values);
-  
-  const unsigned int n_qp_v = qf_v.size();
-  
-
-  typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
-    endc = vspace.get_dh().end();
-  for (; cell!=endc; ++cell) {
-    fe_v.reinit(cell);
-
-    for (unsigned int qp = 0; qp<n_qp_v; ++qp) {
-
-      hardening(cell->index(), qp) = 2e8;
-
-      //plastic_strain(cell->index(), qp)[0][1] = 0;
-
-
-    }
-  }
-
-
-
-}
-
-template <int dim>
-void CamClay<dim>::update_internal_variables(VectorSpace<dim> &vspace,
-                                            Vector<double> &elastic_solution,
-                                             ParsedSymmetricTensorFunction<4, dim> &C)
-{
-
-  //first thing, we need to loop over the cells 
-  //and then the quadrature points
-  QGauss<dim>   qf_v((2*vspace.get_fe().degree) + 1);
-  
-  FEValues<dim> fe_v (vspace.get_fe(), qf_v, 
-                     update_values   | update_gradients |
-                     update_quadrature_points | update_JxW_values);
-  
-  const unsigned int n_qp_v = qf_v.size();
-
-  std::vector<Point<dim> > points(n_qp_v);
-  
-  SymmetricTensor<2,dim> trial_strain;
-  SymmetricTensor<2,dim> trial_stress;
-  SymmetricTensor<2,dim> total_strain;
-
-  std::vector< std::vector<Tensor<1,dim> > > total_grads(n_qp_v, std::vector<Tensor<1,dim> >(dim));
-
-  //average number of NR Iterations per cell 
-  double nr_iterations_ave = 0;
-
-  //This is the material parameter M, for the yield function/surface
-  //right now, I just consider it constant
-  M = 1.0;
-  unsigned int counter = 0;
-  bool marker = false;
-
-  //loop only over the number of cells on this process
-  //we will thread this later
-  typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
-    endc = vspace.get_dh().end();
-  for (; cell!=endc; ++cell) {
-    //reinit with the correct cell
-    fe_v.reinit(cell);
-    int cell_index = cell->index();
-    nr_iterations_ave = 0;
-    marker = false;
-
-    //get the gradients of the PURELY elastic solution
-    fe_v.get_function_gradients(elastic_solution, total_grads);
-
-    //get the quadrature points for the cell
-    points = fe_v.get_quadrature_points();
-
-    //loop over the quadrature points---------------------------------
-    for (unsigned int qp = 0; qp<n_qp_v; ++qp) {
-
-      //get the elastic moduli at this point
-      C_qp = C(points[qp]);
-     
-      //calculate the infinitesimal strain
-      for(unsigned int m=0; m<dim; ++m) {
-       for(unsigned int n=0; n<dim; ++n) {
-         total_strain[m][n] = 0.5*(total_grads[qp][m][n]+total_grads[qp][n][m]);
-       }
-      }
-  
-      //get the trial strain and hardening.  Here we subtract the existing plastic
-      //strain from the new total strain, giving the trial elastic strain
-      //from this step.  The trial hardening is simply the exisiting hardening
-      //parameter
-      trial_strain = (total_strain - plastic_strain(cell_index, qp));
-      //trial_strain *= -1.0;
-
-      //calculate the trial stress.  This is of course just the contraction
-      //of the elasticity tensor and the trial strain
-      trial_stress = C_qp * trial_strain;
-  
-      //if the yield function is greater than zero, solve 
-      //for the new values
-      if ( yield_function(p(trial_stress), q(xi(trial_stress)), 
-                         hardening(cell_index, qp) ) > 0 ) 
-       {
-         marker = true;
-         //deallog << "Cell: " << cell_index
-         //      << " stress01 value: " << trial_stress[0][1] << std::endl; 
-         
-          nr_iterations_ave += solve(cell_index,
-                                    qp,
-                                    total_strain);
-
-       
-       }
-      
-    } //quad points
-    
-    nr_iterations_ave /= n_qp_v;
-    
-    iterations(cell_index) = nr_iterations_ave;
-    
-    if (marker == true) ++counter;
-
-  } //cells
-
-  deallog << "Plastic Deformation in: "
-         << counter << " of " << vspace.get_tria().n_active_cells() << " cells." << std::endl;
-  
-}
-
-template <int dim>
-inline double CamClay<dim>::p(const SymmetricTensor<2, dim> &stress)
-{
-
-  //expicitly assumes specific 2D symmetry - NO SHEAR 
-  double pval = (2*stress[0][0] + stress[1][1]);
-
-  pval *= -1.0/3.0;
-
-  return pval;
-
-}
-
-template <int dim>
-inline SymmetricTensor<2,3> CamClay<dim>::xi(const SymmetricTensor<2,dim> &stress)
-{
-
-  SymmetricTensor<2,3> xival;
-
-  xival[0][0] = stress[0][0];
-  xival[1][1] = stress[0][0];
-  xival[2][2] = stress[1][1];
-
-  for(unsigned int i=0; i<3; ++i) xival[i][i] += p(stress);
-
-  return xival;
-}
-
-template <int dim>
-inline double CamClay<dim>::q(const SymmetricTensor<2, 3> &xi)
-{
-
-  double qval = std::sqrt(3.0/2.0) * xi.norm();
-
-  return qval;
-
-}
-
-template <int dim>
-double CamClay<dim>::solve(const int &index,
-                           const unsigned int &qp,
-                           const SymmetricTensor<2,dim> &total_strain)
-{
-  /*What we have is a system of non-linear ODE that need to
-    be solved for the proper strains, plastic multiplier, and
-    hardening coefficient.  The plastic multiplier is relevant only to
-    this function, so we will only see it here*/
-  
-  /*This version is simpler - it is NOT using the invariant based
-    methods of Claudio Tamagnini*/
-  
-  //these numbers are the 
-  //stress components+hardening+plastic multiplier
-  unsigned int size =( (dim*dim) + dim + 4 )/2;
-  
-  //we need vectors and a matrix
-  Vector<double> cell_res(size);
-  Vector<double> cell_sol(size);
-  Vector<double> cell_prev(size);
-  FullMatrix<double> cell_jac(size,size);
-  
-  //these will be needed later
-  double first_norm = 0;
-  unsigned int iters = 0;
-  
-  //create the guess, we start with the previously converged values of
-  //the plastic strain and hardening, and zero for the plastic mult
-  if (dim == 1) {
-    cell_sol(0) = plastic_strain(index, qp)[0][0]; //e_p
-    cell_sol(1) = hardening(index, qp); //k
-    cell_sol(2) = 0; //delta gamma
-  }
-  if (dim == 2) { 
-    cell_sol(0) = plastic_strain(index, qp)[0][0]; //e_p 11
-    cell_sol(1) = plastic_strain(index, qp)[1][1]; //e_p 22
-    cell_sol(2) = plastic_strain(index, qp)[0][1]; //e_p 12/21
-    cell_sol(3) = hardening(index, qp); //k
-    cell_sol(4) = 0; //delta gamma
-  }
-  if (dim == 3) {
-    cell_sol(0) = plastic_strain(index, qp)[0][0]; //e_p xx
-    cell_sol(1) = plastic_strain(index, qp)[1][1]; //e_p yy
-    cell_sol(2) = plastic_strain(index, qp)[2][2]; //e_p zz
-    cell_sol(3) = plastic_strain(index, qp)[1][2]; //e_p yz
-    cell_sol(4) = plastic_strain(index, qp)[0][2]; //e_p xz
-    cell_sol(5) = plastic_strain(index, qp)[0][1]; //e_p xy
-    cell_sol(6) = hardening(index, qp); //k
-    cell_sol(7) = 0; //delta gamma 
-  }
-  //since the guess of the solution is the previously converged values
-  cell_prev = cell_sol;
-  first_norm = 0;
-  //loop over a newton raphson scheme.
-  //for right now, max nr iterations is set to 10
-  for(unsigned int n=0; n<10; ++n) {
-    cell_res = 0;
-    cell_jac = 0; 
-
-    //compute the residual using the guess
-    compute_res_jac(cell_res, cell_jac, 
-                   cell_sol, cell_prev, total_strain); 
-    
-    cell_res *= -1.0;
-
-    if ( (index == 20) && (qp == 0)) {
-      deallog << "Cell Index: " << index
-            << " qp: " << qp << " Iters: " << iters 
-            << " Rnorm: " << cell_res.linfty_norm() << std::endl;
-      
-      // for (unsigned int i=0; i<size; ++i) 
-      //deallog << "     Residual Comp: " << cell_res(i)  << std::endl;
-    }
-
-    
-    //check for convergence, for right now just use 1e-10
-    if(n == 0) first_norm = cell_res.linfty_norm();
-    if((cell_res.linfty_norm()/first_norm) < 1e-13) break;
-    
-    //if we didn't converge, we are doing an iteration
-    ++iters;
-
-
-     //deallog << "Solution: " << cell_sol(size-2) << std::endl;
-
-     // std::ofstream out("matrix.txt");    
-
-    //cell_jac.print_formatted(out,3,true,0,"0",1,0);
-
-    //invert the 5x5 full matrix
-    cell_jac.gauss_jordan();
-    
-    //vmult the inverted matrix and add it to the solution
-    cell_jac.vmult_add(cell_sol, cell_res);
-    
-  }
-
-  //put the values into the proper places
-  for(unsigned int i=0; i<dim; ++i) {
-    for(unsigned int j=0; j<dim; ++j) {
-      if (i >= j) plastic_strain(index, qp)[i][j] = cell_sol(sym2voigt(i,j));
-    }
-  }
-  hardening(index,qp) = cell_sol(size-2);
-  iter_table(index,qp) = iters;  
-
-  return iters;
-}
-
-template <int dim>
-void CamClay<dim>::compute_res_jac(Vector<double> &res,
-                                  FullMatrix<double> &jac,
-                                  const Vector<double> &sol,
-                                  const Vector<double> &prev_sol,
-                                  const SymmetricTensor<2,dim> &total_strain)
-{
-  //the critical function, computing the residual and the jacobian
-  //basically, if we can assemble the residual using the sacado
-  //doubles, then we have done all of the work.  We just need to be careful about
-  //NOT exploiting the fact that the total strain is symmetric.
-
-  typedef Sacado::Fad::DFad<double> fad_double;
-  unsigned int size = sol.size();
-
-  std::vector< std::vector<fad_double> > estrain(dim, std::vector<fad_double>(dim));
-
-  //vectors of the unknowns and independent vars
-  std::vector<fad_double> x(size); //this is ep11,ep22,ep12,k,dgamma
-  std::vector<fad_double> off_diag(((dim*dim)-dim)/2);
-  std::vector<fad_double> R(size); //this is ep11,ep22,ep12,k,dgamma
-
-  for (unsigned int i=0; i<(2*(size-1))-dim; ++i) {
-    //assign the values of the independent variables
-    if( i<size) {
-      x[i] = sol(i);
-      x[i].diff(i,(2*(size-1))-dim);
-    }
-
-    if(i>=size) {
-      off_diag[i-size] = sol(i-3);
-      off_diag[i-size].diff(i,(2*(size-1))-dim);
-    }
-
-  }
-
-  for(unsigned int i=0; i<dim; ++i) {
-    for (unsigned int j=0; j<dim; ++j) {
-
-      if(i >= j) estrain[i][j] =  (total_strain[i][j] - x[sym2voigt(i,j)]);
-
-      if(i < j) estrain[i][j] = (total_strain[i][j] - off_diag[sym2voigt(i,j)]);
-
-    }
-  }
-
-  //now we are set, we have estrain, with the right values, 
-  //but with different sacado variables, and useable notation
-  std::vector< std::vector<fad_double> > stress(dim, std::vector<fad_double>(dim));
-  fad_double p;
-  p=0;
-  for(unsigned int i=0; i<dim; ++i) {
-    for(unsigned int j=0; j<dim; ++j) {
-      for (unsigned int m=0; m<dim; ++m) {
-       for (unsigned int n=0; n<dim; ++n) {
-         
-         stress[i][j] += C_qp[i][j][m][n] * estrain[m][n];
-         
-       }
-      }
-
-      if (i==j) p += stress[i][j];
-
-    }
-  }
-  
-  p *= (-1.0/double(dim));
-
-  //deallog << "P value: " << p.val() <<std::endl;
-
-  std::vector< std::vector<fad_double> > s(dim, std::vector<fad_double>(dim));
-  fad_double q;
-  fad_double trs;
-  trs = 0;
-
-  for (unsigned int i=0; i<dim; ++i) {
-    for (unsigned int j=0; j<dim; ++j) {
-
-      s[i][j] = stress[i][j];
-
-      if (i == j) {
-
-       s[i][j] += p;
-       trs += s[i][j];
-       
-      }
-
-       q += s[i][j]*s[i][j];
-
-    }
-  }
-
-  q = std::sqrt(3.0*q/2.0);
-  //deallog << " s00 value: " << estrain[0][0].val() << std::endl;
-  //deallog << " s01 value: " << estrain[1][1].val() << std::endl;
-  //deallog << " q value: " << q.val() << std::endl;
-  //deallog << " dg value: " << x[size-1].val() << std::endl;
-
-  for(unsigned int i=0; i<(size-2); ++i) {
-
-    R[i] = x[i] - prev_sol(i) - ( x[size-1]*3.0*s[voigt2sym(i)[0]][voigt2sym(i)[1]]/M/M );
-
-    //deallog << R[i] << std::endl;
-
-    if(i<dim) R[i] -= x[size-1]*( ((x[size-2] - 2.0*p)/double(dim)) + (3.0*trs/M/M/double(dim)) ); 
-
-  }
-
-  //if (R[0].fastAccessDx(0) - R[1].fastAccessDx(1) > 1e-8) deallog << "Shit is Broke" <<std::endl;
-
-  R[size-2] = x[size-2] - prev_sol(size-2) - (16.0*x[size-2]*x[size-1]*(2.0*p - x[size-2]));
-
-  R[size-1] = (q*q/M/M) + (p*(p-x[size-2]));
-
-  for (unsigned int i=0; i<size; ++i) {
-
-    res(i) = R[i].val();
-    //deallog << " R value: " << R[i] << std::endl;
-
-    for (unsigned int j=0; j<size; ++j) {
-      
-      jac(i,j) = R[i].fastAccessDx(j);
-
-      
-    }//j
-    
-  }//i
-
-
-}
-
-
-template <int dim>
-inline std::vector<unsigned int> CamClay<dim>::voigt2sym(const unsigned int i)
-{
-
-  std::vector<unsigned int> return_indices(2);
-
-  if (dim == 2) {
-    if (i == 0) { return_indices[0]=0; return_indices[1]=0; }
-    if (i == 1) { return_indices[0]=1; return_indices[1]=1; }
-    if (i == 2) { return_indices[0]=0; return_indices[1]=1; }
-  }
-
-  if (dim == 3) {
-    //finish this at some time
-  }
-
-  return return_indices;
-
-}
-
-template <int dim>
-inline double CamClay<dim>::delta(const unsigned int i,
-                                 const unsigned int j)
-{
-  double delta_val=0;
-
-  if(i == j) delta_val = 1;
-
-  return delta_val;
-
-}
-
-template <int dim>
-inline double CamClay<dim>::delta(const unsigned int i)
-{
-  double delta_val=1;
-
-  if(i >= dim) delta_val = 0;
-
-  return delta_val;
-
-}
-
-
-template <int dim>
-inline unsigned int CamClay<dim>::sym2voigt(const unsigned int i,
-                                           const unsigned int j)
-{
-  unsigned int voigt_value = 0;
-
-  if (dim == 2) {
-
-    if ( (i == 0) && (j == 0) ) voigt_value = 0;
-    if ( (i == 1) && (j == 1) ) voigt_value = 1;
-    if ( (i == 0) && (j == 1) ) voigt_value = 2;
-    if ( (i == 1) && (j == 0) ) voigt_value = 2;
-
-
-  }
-
-  if (dim == 3) {
-    //needs finished
-
-  }
-
-  return voigt_value;
-
-}
-
-template <int dim>
-double CamClay<dim>::h(const double p,
-                       const double k)
-{
-  
-  double hval = k * ( 2*p - k ); 
-  
-  return hval;
-  
-}
-
-template <int dim>
-Vector<double> CamClay<dim>::dF_dstress(const SymmetricTensor<2, dim> &stress,
-                                        const double k)
-{
-  SymmetricTensor<2,dim> tmp;
-  Vector<double> tmpv(((dim*dim)+dim)/2);
-  
-  double p = (1.0/3.0) * first_invariant(stress);
-  
-  tmp = stress;
-  for(unsigned int i=0; i<dim; ++i) tmp[i][i] -= p;
-  
-  double norm_xi = tmp.norm();
-  
-  double q = (std::sqrt(2.0/3.0)) * norm_xi;
-  
-  tmp = stress;
-  
-  tmp *= (q/M/M) * std::sqrt(8.0/3.0) * (1.0/norm_xi);
-  
-  for(unsigned int i=0; i<dim; ++i )
-    tmp[i][i] += (2.0/3.0)*p - (1.0/3.0)*k - (q*p/M/M)*std::sqrt(8.0/3.0)*(1.0/norm_xi);
-  
-  tmpv(0) = tmp[0][0];
-  tmpv(1) = tmp[1][1];
-  tmpv(2) = tmp[0][1];
-  
-  return tmpv;
-  
-}
-
-template <int dim>
-inline double CamClay<dim>::yield_function(const double p, const double q,
-                                          const double k)
-{
-  double f_value = p*(p-k) + (q*q/M/M);
-
-  return f_value;
-
-}
-
-
-template <int dim>
-void CamClay<dim>::project_strain(ParameterHandler &prm,
-                                 VectorSpace<dim> &vspace,
-                                 SparseMatrix<double> &A)
-{
-
-  LocalAssemblePlasticProject<dim> local_plastic_project;
-  local_plastic_project.reinit(vspace.get_fe(), plastic_strain);
-  local_plastic_project.parameters(prm);
-
-  MyTools::assemble_rhs(vspace.get_dh(), vspace.get_hang(),
-                       -1, solution, local_plastic_project);
-
-  SparseDirectUMFPACK  direct_solver;
-
-  direct_solver.initialize(A);
-
-  direct_solver.solve(solution);
-
-  vspace.get_hang().distribute(solution);
-
-}
-
-template <int dim>
-void CamClay<dim>::build_matrix(VectorSpace<dim> &vspace)
-{
-  QGauss<dim> quad((2*vspace.get_fe().degree) + 1);
-
-  MatrixCreator::create_mass_matrix(vspace.get_dh(), quad, MM);
-
-}
-
-template <int dim>
-void CamClay<dim>::project_hardening(VectorSpace<dim> &vspace)
-{
-  //we need a mass matrix for a scalar variable over the mesh
-  //so, lets build one from scratch
-
-  LocalAssembleScalarProject<dim>  local_ass;
-  local_ass.reinit(vspace.get_fe(), hardening, iter_table);
-
-  MyTools::assemble_rhs(vspace.get_dh(), vspace.get_hang(),
-                       -1, sol_hard_iter, local_ass);
-
-  SparseDirectUMFPACK  direct_solver;
-
-  direct_solver.initialize(MM);
-
-  direct_solver.solve(sol_hard_iter);
-
-  vspace.get_hang().distribute(sol_hard_iter);
-
-  //project the iterations while we are at it
-  //local_ass.reinit(vspace_cc, iter_table);
-  //MyTools::assemble_rhs(vspace_cc.get_dh(), vspace_cc.get_hang(),
-  //                   -1, sol_iterations, local_ass);
-
-
-
-}
-
diff --git a/deal.II/examples/step-41/include/domain.h b/deal.II/examples/step-41/include/domain.h
deleted file mode 100644 (file)
index e4967c0..0000000
+++ /dev/null
@@ -1,72 +0,0 @@
-#ifndef DOMAIN_H
-#define DOMAIN_H
-#include <fstream>
-
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/base/path_search.h>
-#include <deal.II/fe/mapping.h>
-#include <deal.II/dofs/dof_handler.h>
-
-using namespace dealii;
-
-/** 
-  Domain object. 
-*/
-template <int dim>
-class Domain : public Subscriptor
-{
-  public:
-      /** Empty constructor. */
-      Domain ();
-
-      /** Full constructor. */
-      Domain (ParameterHandler &prm);
-
-      ~Domain ();
-
-      /** Reinit. */
-      void reinit(ParameterHandler &prm);
-
-      /** Read mesh file name, etc. */
-      void parse_parameters(ParameterHandler &prm);
-
-      /** Generate entries in the given parameter file. */
-      static void declare_parameters(ParameterHandler &prm);
-
-      /** Generate the mesh. In this program, the mesh can be read from an
-       input file generated with gmsh (http://www.geuz.org/gmsh/). */
-      void create_mesh ();
-      
-      /** Write the mesh. */
-      void output_mesh(std::ostream &out) const;
-      
-
-      /** Write the mesh on the file specified by the parameter handler. */
-      void output_mesh() const;
-
-      /** Reference to the triangulation. */
-      inline Triangulation<dim> & get_tria() {
-       Assert(initialized, ExcNotInitialized());
-       return tria;
-      }
-
-  private:
-      bool initialized;
-
-      bool read_mesh;
-
-      PathSearch search_mesh;
-
-      std::string input_mesh_file_name;
-      std::string input_mesh_format;
-      std::string output_mesh_file_name;
-      
-      /** Holds the coarse triangulation. */
-      Triangulation<dim>   tria;
-      /** Helper class to output the grid. */
-      GridOut gridout;
-};
-#endif
diff --git a/deal.II/examples/step-41/include/domain.templates.h b/deal.II/examples/step-41/include/domain.templates.h
deleted file mode 100644 (file)
index 188f006..0000000
+++ /dev/null
@@ -1,120 +0,0 @@
-#include "domain.h"
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_in.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/fe/fe.h>
-
-using namespace std;
-
-template <int dim> 
-Domain<dim>::Domain() :
-    search_mesh("MESH", 1)
-{
-    initialized = false;
-}
-
-template <int dim>
-Domain<dim>::Domain(ParameterHandler &prm) :
-    search_mesh("MESH", 1)
-{
-    reinit(prm);
-}
-
-template <int dim> 
-Domain<dim>::~Domain()
-{
-  //     tria.clear();
-}
-
-
-template <int dim>
-void Domain<dim>::reinit(ParameterHandler &prm) 
-{
-    deallog.push("DOMAIN");
-    parse_parameters(prm);
-    deallog << "Generating coarse triangulation." << endl;
-    create_mesh();
-    output_mesh();
-    deallog.pop();
-}
-
-
-template <int dim>
-void Domain<dim>::declare_parameters(ParameterHandler &prm) 
-{
-    prm.enter_subsection("Domain Parameters");
-    prm.declare_entry ("Read domain mesh from file", "false", Patterns::Bool(),
-                      "If this is false, then the input mesh file below is ignored and a hyper-cube is created.");
-    prm.declare_entry ("Path of domain mesh files", "mesh/", Patterns::Anything());
-    prm.declare_entry ("Input mesh file", "square", Patterns::Anything()); 
-    prm.declare_entry ("Input mesh format", "ucd", 
-                      Patterns::Selection(GridIn<dim>::get_format_names())); 
-    prm.declare_entry ("Output mesh file", "square_out", Patterns::Anything()); 
-
-    prm.enter_subsection("Grid Out Parameters");
-    GridOut::declare_parameters(prm);
-    prm.leave_subsection();
-
-    prm.leave_subsection();
-
-}
-
-    template <int dim>
-void Domain<dim>::parse_parameters(ParameterHandler &prm) 
-{
-    prm.enter_subsection("Domain Parameters");
-    read_mesh          = prm.get_bool ("Read domain mesh from file");
-    search_mesh.add_path(prm.get ("Path of domain mesh files"));
-    input_mesh_file_name = prm.get ("Input mesh file");
-    output_mesh_file_name = prm.get ("Output mesh file");
-    input_mesh_format  = prm.get ("Input mesh format");
-
-    prm.enter_subsection("Grid Out Parameters");
-    gridout.parse_parameters(prm);
-    prm.leave_subsection();
-
-    prm.leave_subsection();
-}
-
-template <int dim>
-void Domain<dim>::create_mesh() 
-{
-  if(read_mesh) {
-    GridIn<dim> grid_in;
-    grid_in.attach_triangulation (tria);
-    string mfilen = search_mesh.find
-      (input_mesh_file_name, 
-       grid_in.default_suffix(grid_in.parse_format(input_mesh_format)), 
-       "r");
-    ifstream mfile(mfilen.c_str());
-    grid_in.read(mfile, GridIn<dim>::parse_format(input_mesh_format));
-  } else {
-      Point<dim> corner;
-      for(unsigned int d=0; d<dim; ++d) corner[d] = 1.;
-    GridGenerator::hyper_rectangle (tria, Point<dim>(), corner, true);
-//        GridTools::partition_triangulation (n_mpi_processes, tria);
-  }
-  initialized = true;
-}
-  
-
-template <int dim>
-void Domain<dim>::output_mesh(std::ostream &out) const {
-  Assert(initialized, ExcNotInitialized());
-  gridout.write (tria, out);
-}
-
-
-template <int dim>
-void Domain<dim>::output_mesh() const {
-  if(gridout.default_suffix() != "") {
-    std::ofstream out_file ((output_mesh_file_name + gridout.default_suffix()).c_str());
-    output_mesh(out_file);
-    out_file.close();
-  }
-}
diff --git a/deal.II/examples/step-41/include/error_handler.h b/deal.II/examples/step-41/include/error_handler.h
deleted file mode 100644 (file)
index fbc52b6..0000000
+++ /dev/null
@@ -1,177 +0,0 @@
-#ifndef LH_ERROR_HANDLER_H
-#define LH_ERROR_HANDLER_H
-
-#include <fstream>
-
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_constraints.h>
-#include <deal.II/lac/vector.h>
-
-#include <deal.II/grid/tria.h>
-
-// #include <deal.II/numerics/error_estimator.h>
-#include <deal.II/base/function.h>
-#include <deal.II/numerics/solution_transfer.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <deal.II/base/convergence_table.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/parameter_handler.h>
-
-#include <map>
-
-enum NormFlags { 
-    None = 0x00, 
-    Linfty = 0x01,
-    L2 = 0x02, 
-    W1infty = 0x04,
-    H1 = 0x08,
-    AddUp = 0x10
-};
-
-using namespace dealii;
-
-template <int dim, typename VEC=Vector<double> >
-class ErrorHandler : public Subscriptor
-{
-public:
-    /** The constructor takes the mpi initialization stuff. */
-    ErrorHandler ();
-    
-    /** Initialize the given values for the paramter file. */
-    static void declare_parameters(ParameterHandler &prm, 
-                                   unsigned int ntables=1);
-    
-    /** Parse the given parameter handler. */
-    void parse_parameters(ParameterHandler &prm);
-  
-    /** Calculate the error of the numeric solution in variuous norms. Store
-        the result in the given table. */
-    void error_from_exact(const DoFHandler<dim> & vspace, 
-                          const VEC &solution, 
-                          const Function<dim> &exact,
-                          unsigned int table_no = 0,
-                          double dt=0.); 
-  
-    /** Difference between two solutions in two different vector spaces. */
-    void difference(const DoFHandler<dim> &, const VEC &,
-                    const DoFHandler<dim> &, const VEC &, 
-                    unsigned int table_no = 0, double dt=0.); 
-  
-    /** Difference between two solutions in the same vector space. */
-    void difference(const DoFHandler<dim> &, const VEC &, 
-                    const VEC &, unsigned int table_no = 0, double dt=0.);
-    
-    /** By default output first table. */
-    void output_table(const unsigned int table_no=0);
-    
-private:  
-    /** Error results.*/
-    std::vector<ConvergenceTable>  tables;  
-    
-    /** Headers for tables and output. Contains the name of the solution
-        components. */
-    std::vector<std::string> headers;
-    
-    /** Headers for latex tables. Contains the name of the solution
-        components. */
-    std::vector<std::string> latex_headers;
-    
-    /** Captions for latex. */
-    std::vector<std::string> latex_captions;
-    
-    /** Names of the tables. */
-    std::vector<std::string> names;
-    
-    /** Type of error to compute per components. */
-    std::vector<std::vector<NormFlags> > types;
-    
-    /** The parameters have been read. */
-    bool initialized;
-    
-    /** Compute the error. If this is false, all functions regarding
-        errors are disabled and don't do anything.*/
-    bool compute_error;
-    
-    /** Add convergence rates. */
-    std::vector<bool> add_rates;
-    
-    /** Write the error files. */
-    bool write_error;
-    
-    /** Output the error file also on screen. */
-    bool output_error;
-    
-    /** The error file format. */
-    std::string error_file_format;
-    
-    /** The extra column to add to the tables. */    
-    std::vector<std::map<std::string, bool> > extras;
-
-    /** Wether or not to calculate the rates according to the given keys. */
-    std::vector<std::string> rate_keys;
-};
-
-/**
- * Global operator which returns an object in which all bits are set
- * which are either set in the first or the second argument. This
- * operator exists since if it did not then the result of the bit-or
- * <tt>operator |</tt> would be an integer which would in turn trigger
- * a compiler warning when we tried to assign it to an object of type
- * NormFlags.
- */
-inline
-NormFlags
-operator | (NormFlags f1, NormFlags f2)
-{
-  return static_cast<NormFlags> (
-                                static_cast<unsigned int> (f1) |
-                                static_cast<unsigned int> (f2));
-}
-
-/**
- * Global operator which sets the bits from the second argument also
- * in the first one.
- */
-inline
-NormFlags &
-operator |= (NormFlags &f1, NormFlags f2)
-{
-  f1 = f1 | f2;
-  return f1;
-}
-
-
-/**
- * Global operator which returns an object in which all bits are set
- * which are set in the first as well as the second argument. This
- * operator exists since if it did not then the result of the bit-and
- * <tt>operator &</tt> would be an integer which would in turn trigger
- * a compiler warning when we tried to assign it to an object of type
- * NormFlags.
- */
-inline
-NormFlags
-operator & (NormFlags f1, NormFlags f2)
-{
-  return static_cast<NormFlags> (
-                                static_cast<unsigned int> (f1) &
-                                static_cast<unsigned int> (f2));
-}
-
-
-/**
- * Global operator which clears all the bits in the first argument if
- * they are not also set in the second argument.
- */
-inline
-NormFlags &
-operator &= (NormFlags &f1, NormFlags f2)
-{
-  f1 = f1 & f2;
-  return f1;
-}
-
-#endif
diff --git a/deal.II/examples/step-41/include/error_handler.templates.h b/deal.II/examples/step-41/include/error_handler.templates.h
deleted file mode 100644 (file)
index 2b646b6..0000000
+++ /dev/null
@@ -1,396 +0,0 @@
-#include "../include/error_handler.h"
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/grid/grid_tools.h>
-
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/fe/fe.h>
-#include <cstdio>
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include <string>
-
-template <int dim, typename VECTOR>
-ErrorHandler<dim,VECTOR>::ErrorHandler ()
-{
-    initialized = false;
-}
-
-template <int dim, typename VECTOR>
-void ErrorHandler<dim,VECTOR>::declare_parameters (ParameterHandler &prm, 
-                                                   unsigned int ntables)
-{
-    prm.declare_entry ("Write error files", "false", Patterns::Bool());
-    prm.declare_entry ("Output error tables", "true", Patterns::Bool());
-    prm.declare_entry ("Error file format", "tex", Patterns::Selection("tex|txt"));  
-    prm.declare_entry ("Compute error", "true", Patterns::Bool());
-    prm.declare_entry ("Table names", "error", Patterns::Anything(),
-                       "Comma separated list of table names. ");
-    prm.declare_entry ("Solution names", "u", Patterns::Anything(), 
-                       "Comma separated list of names for the components. This "
-                       "will be used both for error tables in text format and to "
-                       "output the solution to a file. Note that in the case "
-                       "of a vector function the error name which is used to "
-                       "compute the norm (supposing the type of the other "
-                       "components is 'Add') is the first one.");
-    prm.declare_entry ("Solution names for latex", "u", Patterns::Anything(),
-                       "Comma separated version of the same thing as above for "
-                       "the latex version of the table.");
-    
-    // prm.declare_entry ("Ib output format", "msh", Patterns::Selection("raw|msh"));
-    // prm.declare_entry ("Ib input file prefix", "ellipse", Patterns::Anything());
-    for(unsigned int i=0; i<ntables; ++i) {
-        char tmp[10];
-        sprintf(tmp, "Table %d", i);
-        prm.enter_subsection(tmp);
-        
-        prm.declare_entry("List of error norms to compute", "Linfty, L2, H1",
-                          Patterns::Anything(), "Each component is separated by a semicolon, "
-                          "and each norm by a comma. Implemented norms are Linfty, L2, "
-                          "H1 and AddUp, which means that the norm is added to the previous "
-                          "component. Useful for vector valued functions.");
-        prm.declare_entry("Add convergence rates", "true", Patterns::Bool(),
-                          "Evaluate convergence rates and add a column to the table for each "
-                          "computed norm. ");
-        prm.declare_entry("Latex table caption", "error", Patterns::Anything(),
-                          "The caption that will go under the table if we write the file in "
-                          "tex format. The default value for this object is the same name "
-                          "as the table name.");
-        prm.declare_entry("Extra terms", "cells,dofs", 
-                          Patterns::Anything(),
-                          "The extra columns to add to the table.");
-        prm.declare_entry("Rate key", "", 
-                          Patterns::Selection("dofs|cells|dt|"),
-                          "The key to use to compute the convergence rates.");
-        prm.leave_subsection();
-    }
-}
-
-template <int dim, typename VECTOR>
-void ErrorHandler<dim,VECTOR>::parse_parameters (ParameterHandler &prm)
-{
-    write_error = prm.get_bool ("Write error files");
-    output_error = prm.get_bool ("Output error tables");
-    
-    error_file_format = prm.get ("Error file format");  
-    compute_error = prm.get_bool ("Compute error");
-    std::string all_names = prm.get ("Table names");
-    headers = Utilities::split_string_list(prm.get ("Solution names"));
-    latex_headers = Utilities::split_string_list(prm.get ("Solution names for latex"));
-       
-    if (all_names != "") {
-        names = Utilities::split_string_list(all_names);
-        types.resize(names.size(), std::vector<NormFlags> (headers.size()));
-        add_rates.resize(names.size());
-        tables.resize(names.size());
-        latex_captions.resize(names.size());
-        std::map<std::string, bool> extra;
-        extra["dof"] = false;
-        extra["cells"] = false;
-        extra["dt"] = false;
-       
-        extras.resize(names.size(), extra);
-       rate_keys.resize(names.size(), "");
-        
-        for(unsigned int i=0; i<names.size(); ++i) {
-            char tmp[10];
-            sprintf(tmp, "Table %d", i);
-            prm.enter_subsection(tmp);
-            
-            all_names = prm.get("List of error norms to compute");
-            add_rates[i] = prm.get_bool("Add convergence rates");
-           rate_keys[i] = prm.get("Rate key");
-            latex_captions[i] = prm.get("Latex table caption");
-            std::vector<std::string> all_extras = 
-                Utilities::split_string_list(prm.get("Extra terms"));
-            
-            for(unsigned int x=0; x< all_extras.size(); ++x)
-                extras[i][all_extras[x]] = true;
-           
-            prm.leave_subsection();
-            
-            std::vector<std::string> all_comps = Utilities::split_string_list(all_names, ';');
-            // Check that the input string has all the needed fields
-            AssertThrow(all_comps.size() == headers.size(),
-                        ExcDimensionMismatch(all_comps.size() , headers.size()));
-            
-            for(unsigned int j=0; j<all_comps.size(); ++j) {
-                std::vector<std::string> all_types = 
-                    Utilities::split_string_list(all_comps[j]);
-                for(unsigned int k=0; k<all_types.size(); ++k) {
-                    if(all_types[k] == "Linfty") {
-                        types[i][j] |= Linfty;
-                    } else if(all_types[k] == "L2") {
-                        types[i][j] |= L2;
-                    } else if(all_types[k] == "W1infty") {
-                        types[i][j] |= W1infty;
-                    } else if(all_types[k] == "H1") {
-                        types[i][j] |= H1;
-                    } else if(all_types[k] == "AddUp") {
-                        types[i][j] |= AddUp;
-                    } else {
-                        AssertThrow(false, ExcMessage("Didn't recognize a norm type."));
-                    }
-                }
-            }
-        }
-    }
-    initialized = true;
-}
-
-template <int dim, typename VECTOR>
-void ErrorHandler<dim,VECTOR>::output_table (const unsigned int table_no) {
-    if (compute_error) {
-        AssertThrow(initialized, ExcNotInitialized());
-        AssertThrow(table_no < names.size(), ExcIndexRange(table_no, 0, names.size()));
-
-        // Add convergence rates
-        if(add_rates[table_no]) {
-            if(extras[table_no]["dofs"])
-                tables[table_no].omit_column_from_convergence_rate_evaluation("dofs"); 
-            if(extras[table_no]["cells"])
-                tables[table_no].omit_column_from_convergence_rate_evaluation("cells");
-            if(extras[table_no]["dt"])
-                tables[table_no].omit_column_from_convergence_rate_evaluation("dt");
-           if(rate_keys[table_no] == "")
-               tables[table_no].evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
-           else 
-               tables[table_no].evaluate_all_convergence_rates(rate_keys[table_no], ConvergenceTable::reduction_rate_log2);
-        }
-
-        if(output_error) tables[table_no].write_text(std::cout);
-
-        if(write_error) {
-            std::string filename = names[table_no] +
-                "." + error_file_format;
-
-            std::ofstream table_file(filename.c_str());
-
-            if(error_file_format != "txt")
-                tables[table_no].write_tex(table_file);
-            else 
-                tables[table_no].write_text(table_file);
-            table_file.close();
-        }
-    }
-}
-
-template <int dim, typename VECTOR>
-void ErrorHandler<dim,VECTOR>::difference(const DoFHandler<dim> & dh, 
-                                          const VECTOR &solution1,
-                                          const VECTOR &solution2, 
-                                          unsigned int table_no,
-                                          double dt) {
-    AssertThrow(solution1.size() == solution2.size(), ExcDimensionMismatch(
-                    solution1.size(), solution2.size()));
-    VECTOR solution(solution1);
-    solution -= solution2;
-    error_from_exact(dh, solution, 
-                     ConstantFunction<dim>(0, headers.size()), table_no, dt);
-}
-
-
-template <int dim, typename VECTOR>
-void ErrorHandler<dim,VECTOR>::error_from_exact(const DoFHandler<dim> & dh, 
-                                                const VECTOR &solution, 
-                                                const Function<dim> &exact,
-                                                unsigned int table_no,
-                                                double dt) 
-{
-    if (compute_error) {
-        AssertThrow(initialized, ExcNotInitialized());
-        AssertThrow(table_no < types.size(), ExcIndexRange(table_no, 0, names.size()));
-        AssertThrow(exact.n_components == types[table_no].size(), 
-                    ExcDimensionMismatch(exact.n_components, types[table_no].size()));
-
-        deallog.push("Error");
-        deallog << "Calculating Errors." << std::endl;
-        std::vector< std::vector<double> > error( exact.n_components, std::vector<double>(4));
-        const unsigned int n_active_cells = dh.get_tria().n_active_cells();
-        const unsigned int n_dofs=dh.n_dofs();
-
-        if(extras[table_no]["cells"]) {
-            tables[table_no].add_value("cells", n_active_cells);
-            tables[table_no].set_tex_caption("cells", "\\# cells");
-            tables[table_no].set_tex_format("cells", "r");
-        }
-        if(extras[table_no]["dofs"]) {
-            tables[table_no].add_value("dofs", n_dofs);
-            tables[table_no].set_tex_caption("dofs", "\\# dofs");
-            tables[table_no].set_tex_format("dofs", "r");
-        }
-        if(extras[table_no]["dt"]) {
-            tables[table_no].add_value("dt", dt);
-            tables[table_no].set_tex_caption("dt", "\\Delta t");
-            tables[table_no].set_tex_format("dt", "r");
-        }
-
-        bool compute_Linfty = false;
-        bool compute_L2 = false;
-        bool compute_W1infty = false;
-        bool compute_H1 = false;
-        bool add_this = false;
-
-        unsigned int last_non_add = 0;
-
-        for(unsigned int component=0; component < exact.n_components; ++component) {
-            NormFlags norm = types[table_no][component];
-
-            deallog << "Error flags: " << norm << std::endl;
-
-            // Select one Component
-            ComponentSelectFunction<dim> select_component ( component, 1. , exact.n_components);
-
-            Vector<float> difference_per_cell (dh.get_tria().n_active_cells());
-
-            QGauss<dim> q_gauss((dh.get_fe().degree+1) * 2);
-
-            // The add bit is set 
-            add_this = (norm & AddUp);
-
-            if(!add_this) {
-                last_non_add   = component;
-                compute_L2     = ( norm & L2 );
-                compute_H1     = ( norm & H1 );
-                compute_W1infty = ( norm & W1infty ) ;
-                compute_Linfty = ( norm & Linfty );
-            }
-            // if add is set, we do not modify the previous selection
-
-            if(compute_L2) {
-                VectorTools::integrate_difference (//mapping,
-                    dh, //dof_handler,
-                    solution, 
-                    exact,
-                    difference_per_cell,
-                    q_gauss,
-                    VectorTools::L2_norm,
-                    &select_component );
-            }
-
-            const double L2_error = difference_per_cell.l2_norm();
-            difference_per_cell = 0;
-          
-            if(compute_H1) { 
-                VectorTools::integrate_difference (//mapping,
-                    dh, //dof_handler,
-                    solution,
-                    exact,
-                    difference_per_cell,
-                    q_gauss,
-                    VectorTools::H1_norm,
-                    &select_component );
-            }
-            const double H1_error = difference_per_cell.l2_norm();
-            difference_per_cell = 0;
-
-            if(compute_W1infty) {
-                VectorTools::integrate_difference (//mapping,
-                    dh, //dof_handler,
-                    solution,
-                    exact,
-                    difference_per_cell,
-                    q_gauss,
-                    VectorTools::W1infty_norm,
-                    &select_component );
-            }
-           
-            const double W1inf_error = difference_per_cell.linfty_norm();
-
-            if(compute_Linfty) {
-                VectorTools::integrate_difference (//mapping,
-                    dh, //dof_handler,
-                    solution,
-                    exact,
-                    difference_per_cell,
-                    q_gauss,
-                    VectorTools::Linfty_norm,
-                    &select_component );
-            }
-
-            const double Linf_error = difference_per_cell.linfty_norm();
-
-            if(add_this) {
-                AssertThrow(component, ExcMessage("Cannot add on first component!"));
-
-                error[last_non_add][0] = std::max(error[last_non_add][0], Linf_error);
-                error[last_non_add][1] += L2_error;
-                error[last_non_add][2] = std::max(error[last_non_add][2], W1inf_error);
-                error[last_non_add][3] += H1_error; 
-
-            } else {
-
-                error[component][0] = Linf_error;
-                error[component][1] = L2_error;
-                error[component][2] = W1inf_error;
-                error[component][3] = H1_error;
-
-            }
-        }  
-
-        for(unsigned int j=0; j<exact.n_components; ++j) {
-            NormFlags norm = types[table_no][j];
-            // If this was added, don't do anything
-            if( !(norm & AddUp) ) {
-                if(norm & Linfty) {
-                    std::string name = headers[j] + "_Linfty";
-                    std::string latex_name = "$\\| " + 
-                        latex_headers[j] + " - " +
-                        latex_headers[j] +"_h \\|_\\infty $";
-                    double this_error =  error[j][0];
-
-                    tables[table_no].add_value(name, this_error);
-                    tables[table_no].set_precision(name, 3);
-                    tables[table_no].set_scientific(name, true);
-                    tables[table_no].set_tex_caption(name, latex_name);
-                }
-
-                if(norm & L2) {
-                    std::string name = headers[j] + "_L2";
-                    std::string latex_name = "$\\| " + 
-                        latex_headers[j] + " - " +
-                        latex_headers[j] +"_h \\|_0 $";
-                    double this_error =  error[j][1];
-
-                    tables[table_no].add_value(name, this_error);
-                    tables[table_no].set_precision(name, 3);
-                    tables[table_no].set_scientific(name, true);
-                    tables[table_no].set_tex_caption(name, latex_name);
-                }
-                if(norm & W1infty) {
-                    std::string name = headers[j] + "_W1infty";
-                    std::string latex_name = "$\\| " + 
-                        latex_headers[j] + " - " +
-                        latex_headers[j] +"_h \\|_{1,\\infty} $";
-                    double this_error =  error[j][2];
-
-                    tables[table_no].add_value(name, this_error);
-                    tables[table_no].set_precision(name, 3);
-                    tables[table_no].set_scientific(name, true);
-                    tables[table_no].set_tex_caption(name, latex_name);
-                }
-                if(norm & H1){
-                    std::string name = headers[j] + "_H1";
-                    std::string latex_name = "$\\| " + 
-                        latex_headers[j] + " - " +
-                        latex_headers[j] +"_h \\|_1 $";
-                    double this_error =  error[j][3];
-
-                    tables[table_no].add_value(name, this_error);
-                    tables[table_no].set_precision(name, 3);
-                    tables[table_no].set_scientific(name, true);
-                    tables[table_no].set_tex_caption(name, latex_name);
-                }
-            }
-        }
-        deallog.pop();
-    }
-}
-
diff --git a/deal.II/examples/step-41/include/linear_elastic.h b/deal.II/examples/step-41/include/linear_elastic.h
deleted file mode 100644 (file)
index a1e4fe9..0000000
+++ /dev/null
@@ -1,77 +0,0 @@
-#ifndef LINEARELASTIC
-#define LINEARELASTIC
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/base/parsed_function.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/sparsity_pattern.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/sparse_ilu.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/sparse_direct.h>
-#include <deal.II/dofs/dof_tools.h>
-
-
-#include "parsed_symmetric_tensor_function.h"
-#include "local_assemble_elastic_matrix.h"
-#include "local_assemble_elastic_rhs.h"
-#include "vector_space.h"
-#include "my_tools.h"
-
-using namespace dealii;
-using namespace dealii::Functions;
-
-template <int dim>
-class LinearElastic
-{
- public:
-  
-  LinearElastic(); 
-  
-  ~LinearElastic(); 
-
-  void declare_parameters(ParameterHandler &prm);
-
-  void parse_parameters(ParameterHandler &prm);
-
-  void reinit(VectorSpace<dim> &vspace);
-
-  void build_matrix(ParameterHandler &prm,
-                   VectorSpace<dim> &vspace);
-
-  void reinit_step(double &time);
-
-  void build_rhs(VectorSpace<dim> &vspace);
-
-  void solve(VectorSpace<dim> &vspace,
-            double &tolerance);
-
-  SparseMatrix<double> A;
-
-  ParsedSymmetricTensorFunction<4, dim>  C;
-
-  //The solution vectors
-  Vector<double> sol_total;
-  Vector<double> sol_increment;
-                       
- private:
-
-  //various functions for the laplace equation
-  ParsedFunction<dim> dbc;
-  ParsedFunction<dim> nbc;
-  ParsedFunction<dim> bf;
-
-  //The system matrix for linear elasticity
-  SparsityPattern  sp_A;
-  
-  //The distributed right hand side      
-  Vector<double>               rhs;
-
-
-
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/linear_elastic.templates.h b/deal.II/examples/step-41/include/linear_elastic.templates.h
deleted file mode 100644 (file)
index 02b51cf..0000000
+++ /dev/null
@@ -1,160 +0,0 @@
-#include "../include/linear_elastic.h"
-
-template <int dim>
-LinearElastic<dim>::LinearElastic()
-  :
-  dbc(dim),
-  nbc(dim),
-  bf(dim)
-{}
-
-
-template <int dim>
-LinearElastic<dim>::~LinearElastic()
-{
-  A.clear();
-}
-
-template <int dim>
-void LinearElastic<dim>::declare_parameters(ParameterHandler &prm)
-{
-  prm.enter_subsection("Dirichlet Data");
-  Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
-  prm.leave_subsection();
-
-  prm.enter_subsection("f - Body Source");
-  Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
-  prm.leave_subsection();
-
-  prm.enter_subsection("Elastic Moduli");
-  ParsedSymmetricTensorFunction<4, dim>::declare_parameters(prm);
-  prm.leave_subsection();
-
-  prm.enter_subsection("Neumann Data");
-  Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
-  prm.leave_subsection();
-
-}
-
-template <int dim>
-void LinearElastic<dim>::parse_parameters(ParameterHandler &prm)
-{ 
-  prm.enter_subsection("f - Body Source");
-  bf.parse_parameters(prm);
-  prm.leave_subsection();
-
-  prm.enter_subsection("Neumann Data");
-  nbc.parse_parameters(prm);
-  prm.leave_subsection();
-
-  prm.enter_subsection("Elastic Moduli");
-  C.parse_parameters(prm);
-  prm.leave_subsection();
-
-  prm.enter_subsection("Dirichlet Data");
-  dbc.parse_parameters(prm);
-  prm.leave_subsection();
-
-}
-
-template <int dim>
-void LinearElastic<dim>::reinit(VectorSpace<dim> &vspace)
-{
-
-  sp_A.reinit(vspace.n_dofs(),
-             vspace.n_dofs(),
-             vspace.get_dh().max_couplings_between_dofs());
-  
-  DoFTools::make_sparsity_pattern (static_cast<DoFHandler<dim> &> (vspace.get_dh()), 
-                                  sp_A);
-  vspace.get_hang().condense (sp_A);
-  sp_A.compress();
-  
-  A.reinit(sp_A);
-
-  sol_total.reinit(vspace.n_dofs());
-  sol_increment.reinit(vspace.n_dofs());
-  rhs.reinit(vspace.n_dofs());
-
-}
-
-template <int dim>
-void LinearElastic<dim>::build_matrix(ParameterHandler &prm,
-                                     VectorSpace<dim> &vspace)
-{
-
-  LocalAssembleElasticMatrix<dim> local_elastic_matrix;
-  local_elastic_matrix.reinit(vspace.get_fe());
-  local_elastic_matrix.parameters(prm);
-
-  MyTools::assemble(vspace.get_dh(), vspace.get_hang(),
-                   -1, A, local_elastic_matrix);
-
-}
-
-template <int dim>
-void LinearElastic<dim>::reinit_step(double &time)
-{
-  
-  rhs = 0;
-  //sol_increment = 0;
-  bf.set_time(time);
-  nbc.set_time(time);
-  dbc.set_time(time);
-}
-
-template <int dim>
-void LinearElastic<dim>::build_rhs(VectorSpace<dim> &vspace)
-{
-
-  LocalAssembleElasticRHS<dim> elastic_local_rhs;
-
-  elastic_local_rhs.reinit(vspace.get_fe(),
-                          bf, nbc,
-                          vspace.neumann_bc);
-  
-  MyTools::assemble_rhs(vspace.get_dh(), vspace.get_hang(),
-                       -1, rhs, elastic_local_rhs);
-  
-}
-
-template <int dim>
-void LinearElastic<dim>::solve(VectorSpace<dim> &vspace,
-                              double &tolerance)
-{
-
-  //get the DBC map
-  std::map<unsigned int, double>  bv;
-  vspace.interpolate_dirichlet_bc(dbc, bv);
-  
-  //we want to leave the matrix untouched for future projections
-  //so we use a filtered matrix object
-  FilteredMatrix<Vector<double> > filtered_A(A);
-  filtered_A.add_constraints(bv);
-  filtered_A.apply_constraints(rhs, true);
-
-  deallog << "RHS L2 Norm After DBC: " << rhs.l2_norm() << std::endl;
-
-  //make the preconditioner
-  PreconditionJacobi<SparseMatrix<double> > precon;
-  precon.initialize(A, 0.8);
-  FilteredMatrix<Vector<double> > filtered_precon(precon);
-
-  SolverControl control (vspace.n_dofs(),
-                        tolerance*rhs.linfty_norm(),
-                        false, true);
-
-  GrowingVectorMemory<Vector<double> > mem;
-  SolverCG<Vector<double> > solver (control, mem);
-
-  //SparseDirectUMFPACK direct_solver;
-  //direct_solver.initialize(filtered_A);
-
-  solver.solve(filtered_A,sol_total,rhs,precon);
-
-  vspace.get_hang().distribute(sol_total);
-
-  //sol_total += sol_increment;
-
-}
-
diff --git a/deal.II/examples/step-41/include/local_assemble.h b/deal.II/examples/step-41/include/local_assemble.h
deleted file mode 100644 (file)
index 89c0bdb..0000000
+++ /dev/null
@@ -1,66 +0,0 @@
-#ifndef LOCAL_ASSEMBLE_STANDARD
-#define LOCAL_ASSEMBLE_STANDARD
-
-#include "local_assemble_base.h"
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/base/parameter_handler.h>
-
-
-template <int dim, typename DH=MGDoFHandler<dim> >
-  class LocalAssemble : public LocalAssembleBase<dim, DH>
-{
-public:
-#if deal_II_dimension != 1
-  LocalAssemble();
-  
-  /** Given two face finite elements, assemble the two matrices.*/
-  virtual void  assemble_face_terms (FEFaceValuesBase<dim> &fe_v,
-                                    FEFaceValuesBase<dim> &fe_n_v,
-                                    FullMatrix<double> &,
-                                    FullMatrix<double> &);
-
-  virtual void  assemble_face_term
-  (const typename DH::active_cell_iterator& , 
-   const unsigned int,
-   const typename DH::active_cell_iterator& , 
-   const unsigned int,
-   FullMatrix<double> &,
-   FullMatrix<double> &);
-
-  virtual void  assemble_face_term
-  (const typename DH::active_cell_iterator& , 
-   const unsigned int,
-   const unsigned int,
-   const typename DH::active_cell_iterator& , 
-   const unsigned int,
-   FullMatrix<double> &,
-   FullMatrix<double> &);
-
-  virtual void  assemble_face_term
-  (const typename DH::active_cell_iterator& , 
-   const unsigned int,
-   const typename DH::active_cell_iterator& , 
-   const unsigned int,
-   const unsigned int,
-   FullMatrix<double> &,
-   FullMatrix<double> &);
-
-  /** A pointer to fe_face_values objects. */
-  SmartPointer<FEFaceValues<dim> > fe_face_v;
-
-  /** A pointer to fe_face_values objects. */
-  SmartPointer<FEFaceValues<dim> > fe_face_n_v;
-
-  /** A pointer to fe_face_values objects. */
-  SmartPointer<FESubfaceValues<dim> > fe_sub_face_v;
-#endif
-
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble.templates.h b/deal.II/examples/step-41/include/local_assemble.templates.h
deleted file mode 100644 (file)
index 1a2d4a5..0000000
+++ /dev/null
@@ -1,67 +0,0 @@
-#include "local_assemble.h"
-#include <ostream>
-
-#if deal_II_dimension != 1
-
-template <int dim, typename DH>
-LocalAssemble<dim,DH>::LocalAssemble() : 
-  fe_face_v(0, "Local Assemble FeFaceValues Pointer"),
-  fe_face_n_v(0, "Local Assemble FeFaceValues on Neighbor Pointer"),
-  fe_sub_face_v(0, "Local Assemble FesubfaceValues Pointer")
-{
-}
-
-template <int dim, typename DH>
-void  LocalAssemble<dim,DH>::assemble_face_terms (FEFaceValuesBase<dim> &,
-                                                 FEFaceValuesBase<dim> &,
-                                                 FullMatrix<double> &,
-                                                 FullMatrix<double> &)
-{
-  Assert(false, ExcPureFunctionCalled());
-}
-
-  
-template <int dim, typename DH>
-void LocalAssemble<dim,DH>::assemble_face_term
-(const typename DH::active_cell_iterator& cell, 
- const unsigned int face_no,
- const typename DH::active_cell_iterator& n_cell, 
- const unsigned int n_face_no,
- FullMatrix<double> &cell_m,
- FullMatrix<double> &cell_n_m)
-{
-  this->fe_face_v->reinit(cell, face_no);
-  this->fe_face_n_v->reinit(n_cell, n_face_no);
-  assemble_face_terms(*fe_face_v, *fe_face_n_v, cell_m, cell_n_m);
-}
-
-template <int dim, typename DH>
-void LocalAssemble<dim,DH>::assemble_face_term
-(const typename DH::active_cell_iterator& cell, 
- const unsigned int face_no,
- const typename DH::active_cell_iterator& n_cell, 
- const unsigned int n_face_no,
- const unsigned int n_sub_face_no,
- FullMatrix<double> &cell_m,
- FullMatrix<double> &cell_n_m)
-{
-  fe_face_v->reinit(cell, face_no);
-  fe_sub_face_v->reinit(n_cell, n_face_no, n_sub_face_no);
-  assemble_face_terms(*fe_face_v, *fe_sub_face_v, cell_m, cell_n_m);
-}
-
-template <int dim, typename DH>
-void LocalAssemble<dim,DH>::assemble_face_term
-(const typename DH::active_cell_iterator& cell, 
- const unsigned int face_no,
- const unsigned int subface_no,
- const typename DH::active_cell_iterator& n_cell, 
- const unsigned int n_face_no,
- FullMatrix<double> &cell_m,
- FullMatrix<double> &cell_n_m)
-{
-  fe_sub_face_v->reinit(cell, face_no, subface_no);
-  fe_face_v->reinit(n_cell, n_face_no);
-  assemble_face_terms(*fe_sub_face_v, *fe_face_v, cell_m, cell_n_m);
-}
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_base.h b/deal.II/examples/step-41/include/local_assemble_base.h
deleted file mode 100644 (file)
index d8576f3..0000000
+++ /dev/null
@@ -1,102 +0,0 @@
-#ifndef __deal2__local_assemble_base_h
-#define __deal2__local_assemble_base_h
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include "assemble_flags.h"
-
-#include <fstream>
-#include <iostream>
-
-using namespace dealii;
-
-/** Base class for local assemblers. This is the object you should
-  derive your own classes from, if you want to use the VirtualMatrix
-  class or the MyTools::assemble utility. It provides a common
-  interface to the procedure for assembling locally the matrices. All
-  functions are virtual and would throw an exception if called. You
-  should really implement your own class deriving it from this one. An
-  example is provided in the LocalAssembleLaplace class, which
-  assemble the Laplace equations locally for continuous finite
-  elements.
-
-  All methods work basically in the same way. They take a FEValues
-  object and a FullMatrix one as input, and fill the FullMatrix with
-  appropriate values. Exceptions are the face and boundary terms,
-  which take FEFaceValues as input.
-
-  This object is supposed to be used in a Virtual Matrix via the 
-  VirtualMatrix::enter method.
-
-  Here the only thing which is actually done is to create the internal 
-  AssembleFlags object upon construction, which is used to determine which
-  of the methods will be called by the virtual matrix.
- */
-template <int dim, typename DH=MGDoFHandler<dim> >
-class LocalAssembleBase : public Subscriptor
-{
-    public:
-       virtual ~LocalAssembleBase() {};
-       /** This object will be called for each cell of a triangulation. */
-       // virtual void assemble_cell_term(const FEValues<dim>& fe_v,
-       //      FullMatrix<double> &u_v_matrix) const;
-       virtual void assemble_cell_term
-           (const typename DH::active_cell_iterator&, 
-            FullMatrix<double> &cell_m);
-
-       /** This object will be called for each boundary face of a triangulation.*/
-       virtual void assemble_boundary_term
-           (const typename DH::active_cell_iterator&, const unsigned int, 
-            FullMatrix<double> &cell_m);
-
-       /** This object will be called for each face of a
-           triangulation. This one is called when the face is shared
-           by neighbors on the same level of refinement. */
-       virtual void assemble_face_term
-         (const typename DH::active_cell_iterator &cell, 
-          const unsigned int face_no, 
-          const typename DH::active_cell_iterator &neighbor, 
-          const unsigned int neighbor_face_no, 
-          FullMatrix<double> &cell_m,
-          FullMatrix<double> &neighbor_cell_m);
-       
-       /** This object will be called for each face of a
-           triangulation. This one is called when the current face is
-           coarser. */
-       virtual void assemble_face_term
-         (const typename DH::active_cell_iterator &cell, 
-          const unsigned int face_no, 
-          const unsigned int sub_face_no, 
-          const typename DH::active_cell_iterator &neighbor, 
-          const unsigned int neighbor_face_no, 
-          FullMatrix<double> &cell_m,
-          FullMatrix<double> &neighbor_cell_m);
-       
-       /** This object will be called for each face of a
-           triangulation. This one is called when the current face is
-           finer. */
-       virtual void assemble_face_term
-         (const typename DH::active_cell_iterator &cell, 
-          const unsigned int face_no, 
-          const typename DH::active_cell_iterator &neighbor, 
-          const unsigned int neighbor_face_no, 
-          const unsigned int neighbor_subface_no,
-          FullMatrix<double> &cell_m,
-          FullMatrix<double> &neighbor_cell_m);
-       
-       /** Assemble rhs. This will be called for each cell.*/
-       virtual void assemble_rhs_term
-         (const typename DH::active_cell_iterator&, 
-          Vector<double> &);
-
-       /** This object will be called for each boundary face of a triangulation.*/
-       virtual void assemble_rhs_boundary_term
-         (const typename DH::active_cell_iterator&, const unsigned int, 
-          Vector<double> &);
-       
-       AssembleFlags flags;
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_base.templates.h b/deal.II/examples/step-41/include/local_assemble_base.templates.h
deleted file mode 100644 (file)
index 98fdc4e..0000000
+++ /dev/null
@@ -1,68 +0,0 @@
-#include "local_assemble_base.h"
-                
-template <int dim, typename DH>  
-void LocalAssembleBase<dim,DH>::assemble_cell_term (const typename DH::active_cell_iterator&, 
-                                                FullMatrix<double> &)
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-}
-
-template <int dim, typename DH>  
-void LocalAssembleBase<dim,DH>::assemble_boundary_term
-(const typename DH::active_cell_iterator&,const  unsigned int, 
- FullMatrix<double> &) 
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-}
-
-template <int dim, typename DH>  
-void LocalAssembleBase<dim,DH>::assemble_face_term
-(const typename DH::active_cell_iterator&,const unsigned int, 
- const typename DH::active_cell_iterator&,const unsigned int, 
- FullMatrix<double> &,
- FullMatrix<double> &) 
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-}
-
-template <int dim, typename DH>  
-void LocalAssembleBase<dim,DH>::assemble_face_term
-(const typename DH::active_cell_iterator&,
- const unsigned int, const unsigned int,
- const typename DH::active_cell_iterator&,const unsigned int, 
- FullMatrix<double> &,
- FullMatrix<double> &) 
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-}
-
-template <int dim, typename DH>  
-void LocalAssembleBase<dim,DH>::assemble_face_term
-(const typename DH::active_cell_iterator&,const unsigned int, 
- const typename DH::active_cell_iterator&,
- const unsigned int, 
- const unsigned int,
- FullMatrix<double> &,
- FullMatrix<double> &) 
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-}
-
-
-template <int dim, typename DH>  
-void LocalAssembleBase<dim,DH>::assemble_rhs_term
-(const typename DH::active_cell_iterator&, 
- Vector<double> &)
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-}
-
-template <int dim, typename DH>  
-void LocalAssembleBase<dim,DH>::assemble_rhs_boundary_term
-(const typename DH::active_cell_iterator&, 
- const unsigned, 
- Vector<double> &)
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-}
-
diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_matrix.h b/deal.II/examples/step-41/include/local_assemble_elastic_matrix.h
deleted file mode 100644 (file)
index 64c21f2..0000000
+++ /dev/null
@@ -1,50 +0,0 @@
-#ifndef LOCAL_ASSEMBLE_ELASTIC_MATRIX
-#define LOCAL_ASSEMBLE_ELASTIC_MATRIX
-
-#include "local_assemble.h"
-#include "parsed_symmetric_tensor_function.h"
-#include <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/base/parameter_handler.h>
-
-
-template <int dim>
-class LocalAssembleElasticMatrix : public LocalAssembleBase<dim>
-{
- public:
-  
-  LocalAssembleElasticMatrix(); 
-  
-  ~LocalAssembleElasticMatrix(); 
-  
-  virtual void  assemble_cell_term
-    (const typename MGDoFHandler<dim>::active_cell_iterator& cell, 
-     FullMatrix<double> &);
-
-
-  void reinit(FiniteElement<dim>&);
-
-  void parameters(ParameterHandler &prm);
-                       
- private:
-  /** A pointer to the finite element.*/
-  SmartPointer<FiniteElement<dim> > fe;
-
-  /** A pointer to the rhs function.*/
-  ParsedSymmetricTensorFunction<4, dim>  C;
-
-  /** A pointer to fe_values objects. */
-  SmartPointer<FEValues<dim> > fe_v;
-  
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h b/deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h
deleted file mode 100644 (file)
index 76480c4..0000000
+++ /dev/null
@@ -1,85 +0,0 @@
-#include "../include/local_assemble_elastic_matrix.h"
-#include <deal.II/base/quadrature_lib.h>
-
-template <int dim>
-LocalAssembleElasticMatrix<dim>::LocalAssembleElasticMatrix() : 
-fe(0, "Local Assemble Fe Pointer"),
-  fe_v(0, "Local Assemble FeValues Pointer")
-{
-}
-
-
-template <int dim>
-void LocalAssembleElasticMatrix<dim>::reinit (FiniteElement<dim>& myfe) 
-{ 
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  fe = &myfe;
-  QGauss<dim> quadrature(2*fe->degree + 1);
-  QGauss<dim-1> face_quadrature(2*fe->degree + 1);
-  UpdateFlags flags (update_values |
-                    update_gradients |
-                    update_q_points  |
-                    update_JxW_values);
-  fe_v = new FEValues<dim>(*fe, quadrature, flags);
-    
-  this->flags = assemble_cell;
-    
-}
-
-template <int dim>
-LocalAssembleElasticMatrix<dim>::~LocalAssembleElasticMatrix()
-{
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  fe = 0;
-}
-
-template <int dim>
-void LocalAssembleElasticMatrix<dim>::parameters(ParameterHandler &prm)
-{
-  prm.enter_subsection("Elastic Moduli");
-  C.parse_parameters(prm);
-  prm.leave_subsection();
-}
-
-template <int dim>
-void LocalAssembleElasticMatrix<dim>::assemble_cell_term
-(const typename MGDoFHandler<dim>::active_cell_iterator& cell, 
- FullMatrix<double> &cell_m)
-{
-  cell_m = 0;
-  Assert(fe, ExcNotInitialized());
-  
-  fe_v->reinit(cell);
-
-  std::vector<Point<dim> > points(fe_v->n_quadrature_points);
-  points = fe_v->get_quadrature_points();
-  unsigned int comp_i, comp_j;  
-
-  for(unsigned int point =0; point<fe_v->n_quadrature_points; ++point) {
-    SymmetricTensor<4,dim> C_qp = C(points[point]);
-    
-    for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
-      comp_i = fe->system_to_component_index(i).first;
-      for (unsigned int j=0; j<fe_v->dofs_per_cell; ++j) {
-       comp_j = fe->system_to_component_index(j).first;
-       
-       for(unsigned int b=0; b<dim; ++b) {
-         for(unsigned int n=0; n<dim; ++n) {
-           cell_m(i,j) +=  C_qp[comp_i][b][comp_j][n] * 
-             fe_v->shape_grad(j, point)[n] *
-             fe_v->shape_grad(i, point)[b] *
-             fe_v->JxW(point);
-         } //n
-       }//b
-      } //j
-    }//i
-  }//qp
-}
diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_rhs.h b/deal.II/examples/step-41/include/local_assemble_elastic_rhs.h
deleted file mode 100644 (file)
index a31081e..0000000
+++ /dev/null
@@ -1,65 +0,0 @@
-#ifndef LOCAL_ASSEMBLE_ELASTIC_RHS
-#define LOCAL_ASSEMBLE_ELASTIC_RHS
-
-#include "local_assemble.h"
-#include "parsed_symmetric_tensor_function.h"
-#include <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/base/parameter_handler.h>
-
-
-template <int dim>
-class LocalAssembleElasticRHS : public LocalAssembleBase<dim>
-{
- public:
-  
-  LocalAssembleElasticRHS(); 
-  
-  ~LocalAssembleElasticRHS(); 
-  
-  virtual void assemble_rhs_term
-    (const typename MGDoFHandler<dim>::active_cell_iterator&,
-     Vector<double> &);
-  
-  /** This object will be called for each boundary face of a triangulation.*/
-  virtual void assemble_rhs_boundary_term
-    (const typename MGDoFHandler<dim>::active_cell_iterator&, const unsigned int, 
-     Vector<double> &);
-
-
-  void reinit(FiniteElement<dim>&,
-             Function<dim> &,
-             Function<dim> &,
-             std::map<char, std::vector<bool> > &);
-                       
- private:
-  /** A pointer to the finite element.*/
-  SmartPointer<FiniteElement<dim> > fe;
-
-  /** A pointer to the rhs function.*/
-  SmartPointer<Function<dim> > rhs;
-
-  /** A pointer to the Neumann function.*/
-  SmartPointer<Function<dim> > neumann;
-  
-  /** A map of ids and components for neumann bc. */
-  std::map<char, std::vector<bool> > neumann_map;
-
-  /** A pointer to fe_values objects. */
-  SmartPointer<FEValues<dim> > fe_v;
-  
-  /** A pointer to fe_face_values objects. */
-  SmartPointer<FEFaceValues<dim> > fe_face_v;
-  
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h b/deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h
deleted file mode 100644 (file)
index c428b26..0000000
+++ /dev/null
@@ -1,133 +0,0 @@
-#include "../include/local_assemble_elastic_rhs.h"
-#include <deal.II/base/quadrature_lib.h>
-
-template <int dim>
-LocalAssembleElasticRHS<dim>::LocalAssembleElasticRHS() : 
-fe(0, "Local Assemble Fe Pointer"),
-  fe_v(0, "Local Assemble FeValues Pointer"),
-  fe_face_v(0, "Local Assemble FeFaceValues Pointer")
-{
-}
-
-
-template <int dim>
-void LocalAssembleElasticRHS<dim>::reinit (FiniteElement<dim> &myfe,
-                                       Function<dim> &bf,
-                                       Function<dim> &nbc,
-                                       std::map<char, std::vector<bool> > & n_map) 
-{ 
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  if(fe_face_v) {
-    FEFaceValues<dim> * p = fe_face_v;
-    fe_face_v = 0;
-    delete p;
-  }
-  fe = &myfe;
-  QGauss<dim> quadrature(2*fe->degree + 1);
-  QGauss<dim-1> face_quadrature(2*fe->degree + 1);
-  UpdateFlags flags (update_values |
-                    update_gradients |
-                    update_q_points  |
-                    update_JxW_values);
-  fe_v = new FEValues<dim>(*fe, quadrature, flags);
-    
-  fe_face_v = new FEFaceValues<dim>(*fe, face_quadrature, flags);
-    
-  this->flags = 
-    assemble_rhs_cell|
-    assemble_rhs_boundary;
-    
-  rhs = &bf;
-  neumann = &nbc;
-  neumann_map = n_map;
-}
-
-template <int dim>
-LocalAssembleElasticRHS<dim>::~LocalAssembleElasticRHS()
-{
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  if(fe_face_v) {
-    FEFaceValues<dim> * p = fe_face_v;
-    fe_face_v = 0;
-    delete p;
-  }
-  fe = 0;
-}
-
-
-template <int dim>
-void LocalAssembleElasticRHS<dim>::assemble_rhs_term
-(const typename MGDoFHandler<dim>::active_cell_iterator& cell, 
- Vector<double> &cell_rhs)
-{
-  cell_rhs = 0;
-  Assert(fe, ExcNotInitialized());
-
-  fe_v->reinit(cell);
-  
-  unsigned int size = fe->n_components();
-  unsigned int n_q_points = fe_v->n_quadrature_points;
-  
-  std::vector<Vector<double> > load_vector (n_q_points, Vector<double>(size) );
-  
-  /* Evaluate rhs and solution on the quadrature points. */
-  rhs->vector_value_list (fe_v->get_quadrature_points(), load_vector);
-
-  unsigned int comp_i;
-  for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
-    comp_i = fe->system_to_component_index(i).first;
-    for (unsigned int q_point=0; q_point<n_q_points; ++q_point) {
-      cell_rhs(i) +=
-       ( load_vector[q_point](comp_i) *
-         fe_v->shape_value(i,q_point) *
-         fe_v->JxW(q_point) );
-    }
-  }
-}
-
-template<int dim>
-void LocalAssembleElasticRHS<dim>::assemble_rhs_boundary_term
-(const typename MGDoFHandler<dim>::active_cell_iterator& cell, 
- const unsigned int face_no,
- Vector<double> &cell_rhs)
-{
-  cell_rhs = 0;
-  // See if we need to do anything here
-  char id = cell->face(face_no)->boundary_indicator();
-  if(neumann_map.find(id) == neumann_map.end()) 
-    return;
-  
-  std::vector<bool> & filter = neumann_map[id];
-  
-  Assert(fe, ExcNotInitialized());
-  fe_face_v->reinit(cell,face_no);
-  
-  unsigned int size = fe->n_components();
-  unsigned int n_q_points = fe_face_v->n_quadrature_points;
-  
-  /** Vector of boundary values.*/
-  std::vector<Vector<double> > neumann_vector (n_q_points, Vector<double>(size) );
-  
-  /* Evaluate rhs and solution on the quadrature points. */
-  neumann->vector_value_list (fe_face_v->get_quadrature_points(), neumann_vector);
-  
-  unsigned int comp_i;
-  for (unsigned int i=0; i<fe_face_v->dofs_per_cell; ++i) {
-    comp_i = fe->system_to_component_index(i).first;
-    if(filter[comp_i]) 
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point) {
-       cell_rhs(i) +=
-         ( neumann_vector[q_point](comp_i) *
-           fe_face_v->shape_value(i,q_point) *
-           fe_face_v->JxW(q_point) );
-      }
-  }
-}
diff --git a/deal.II/examples/step-41/include/local_assemble_hardening_project.h b/deal.II/examples/step-41/include/local_assemble_hardening_project.h
deleted file mode 100644 (file)
index fc4be24..0000000
+++ /dev/null
@@ -1,53 +0,0 @@
-#ifndef LOCAL_ASSEMBLE_HARDENING_PROJECT
-#define LOCAL_ASSEMBLE_HARDENING_PROJECT
-
-#include "local_assemble.h"
-#include "parsed_symmetric_tensor_function.h"
-#include <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/base/parameter_handler.h>
-
-
-template <int dim>
-class LocalAssemblePlasticProject : public LocalAssembleBase<dim>
-{
- public:
-  
-  LocalAssemblePlasticProject(); 
-  
-  ~LocalAssemblePlasticProject(); 
-  
-  virtual void assemble_rhs_term
-    (const typename MGDoFHandler<dim>::active_cell_iterator&,
-     Vector<double> &);
-  
-  void parameters(ParameterHandler &prm);
-
-
-  void reinit(FiniteElement<dim>&,
-             Table<2, SymmetricTensor<2,dim> > &);
-                       
- private:
-  /** A pointer to the finite element.*/
-  SmartPointer<FiniteElement<dim> > fe;
-
-  /** A pointer to the plastic strain */
-  Table<dim,SymmetricTensor<2,dim> >  plastic_strain;
-
-  /** A pointer to fe_values objects. */
-  SmartPointer<FEValues<dim> > fe_v;
-
-  ParsedSymmetricTensorFunction<4, dim> C;
-  
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_mass.h b/deal.II/examples/step-41/include/local_assemble_mass.h
deleted file mode 100644 (file)
index 9d3e808..0000000
+++ /dev/null
@@ -1,45 +0,0 @@
-#ifndef LOCAL_ASSEMBLE_MASS_MATRIX
-#define LOCAL_ASSEMBLE_MASS_MATRIX
-
-#include "local_assemble.h"
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-template <int dim, typename DH>
-class LocalAssembleMass : 
-public LocalAssembleBase<dim, DH>
-{
- public:
-  
-  LocalAssembleMass(); 
-  
-  ~LocalAssembleMass(); 
-
-  virtual void  assemble_cell_term
-    (const typename DH::active_cell_iterator& cell, 
-     FullMatrix<double> &);
-  
-  virtual void assemble_rhs_term
-    (const typename DH::active_cell_iterator&,
-     Vector<double> &);
-  
-  void reinit(FiniteElement<dim>&,
-             Function<dim> &);
-       
- private:
-  /** A pointer to the finite element.*/
-  SmartPointer<FiniteElement<dim> > fe;
-
-  /** A pointer to the rhs function.*/
-  SmartPointer<Function<dim> > rhs;
-
-  /** A pointer to fe_values objects. */
-  SmartPointer<FEValues<dim> > fe_v;
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_mass.templates.h b/deal.II/examples/step-41/include/local_assemble_mass.templates.h
deleted file mode 100644 (file)
index da02224..0000000
+++ /dev/null
@@ -1,93 +0,0 @@
-#include "../include/local_assemble_mass.h"
-#include "utilities.h"
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-
-template <int dim, typename DH>
-LocalAssembleMass<dim, DH>::LocalAssembleMass() : 
-    fe(0, "Local Assemble Fe Pointer"),
-    fe_v(0, "Local Assemble FeValues Pointer")
-{
-}
-
-
-template <int dim, typename DH>
-void LocalAssembleMass<dim, DH>::reinit (FiniteElement<dim> &myfe,
-                                          Function<dim> &f) 
-{ 
-    smart_delete(fe_v);
-    fe = &myfe;
-    QGauss<dim> quadrature(2*fe->degree + 1);
-    UpdateFlags flags (update_values |
-                      update_q_points  |
-                      update_JxW_values);
-    fe_v = new FEValues<dim>(*fe, quadrature, flags);
-    
-    this->flags = 
-       assemble_cell|
-       assemble_rhs_cell; 
-    
-    rhs = &f;
-}
-
-template <int dim, typename DH>
-LocalAssembleMass<dim, DH>::~LocalAssembleMass()
-{
-    smart_delete(fe_v);
-    fe = 0;
-}
-
-template <int dim, typename DH>
-void LocalAssembleMass<dim, DH>::assemble_cell_term
-(const typename DH::active_cell_iterator& cell, 
- FullMatrix<double> &cell_m)
-{
-  cell_m = 0;
-  Assert(fe, ExcNotInitialized());
-  
-  fe_v->reinit(cell);
-  unsigned int comp_i = 0, comp_j = 0;
-  for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
-      comp_i = fe->system_to_component_index(i).first;
-      for (unsigned int j=0; j<fe_v->dofs_per_cell; ++j) {      
-         comp_j = fe->system_to_component_index(j).first;
-         if(comp_i == comp_j) 
-             for(unsigned int q_point =0; q_point<fe_v->n_quadrature_points; ++q_point) {
-                 cell_m(i,j) +=  ( fe_v->shape_value(j,q_point) *
-                                   fe_v->shape_value(i,q_point) * 
-                                   fe_v->JxW(q_point) );
-             }
-      }
-  }
-}
-
-
-template <int dim, typename DH>
-void LocalAssembleMass<dim, DH>::assemble_rhs_term
-(const typename DH::active_cell_iterator& cell, 
- Vector<double> &cell_rhs)
-{
-  cell_rhs = 0;
-  Assert(fe, ExcNotInitialized());
-
-  fe_v->reinit(cell);
-  
-  unsigned int size = fe->n_components();
-  unsigned int n_q_points = fe_v->n_quadrature_points;
-  
-  std::vector<Vector<double> > load_vector (n_q_points, Vector<double>(size) );
-  
-  /* Evaluate rhs and solution on the quadrature points. */
-  rhs->vector_value_list (fe_v->get_quadrature_points(), load_vector);
-
-  unsigned int comp_i;
-  for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
-    comp_i = fe->system_to_component_index(i).first;
-    for (unsigned int q_point=0; q_point<n_q_points; ++q_point) {
-      cell_rhs(i) +=
-       ( load_vector[q_point](comp_i) *
-         fe_v->shape_value(i,q_point) *
-         fe_v->JxW(q_point) );
-    }
-  }
-}
diff --git a/deal.II/examples/step-41/include/local_assemble_plastic_project.h b/deal.II/examples/step-41/include/local_assemble_plastic_project.h
deleted file mode 100644 (file)
index 5bd80c8..0000000
+++ /dev/null
@@ -1,53 +0,0 @@
-#ifndef LOCAL_ASSEMBLE_PLASTIC_PROJECT
-#define LOCAL_ASSEMBLE_PLASTIC_PROJECT
-
-#include "local_assemble.h"
-#include "parsed_symmetric_tensor_function.h"
-#include <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/base/parameter_handler.h>
-
-
-template <int dim>
-class LocalAssemblePlasticProject : public LocalAssembleBase<dim>
-{
- public:
-  
-  LocalAssemblePlasticProject(); 
-  
-  ~LocalAssemblePlasticProject(); 
-  
-  virtual void assemble_rhs_term
-    (const typename MGDoFHandler<dim>::active_cell_iterator&,
-     Vector<double> &);
-  
-  void parameters(ParameterHandler &prm);
-
-
-  void reinit(FiniteElement<dim>&,
-             Table<2, SymmetricTensor<2,dim> > &);
-                       
- private:
-  /** A pointer to the finite element.*/
-  SmartPointer<FiniteElement<dim> > fe;
-
-  /** A pointer to the plastic strain */
-  Table<dim,SymmetricTensor<2,dim> >  plastic_strain;
-
-  /** A pointer to fe_values objects. */
-  SmartPointer<FEValues<dim> > fe_v;
-
-  ParsedSymmetricTensorFunction<4, dim> C;
-  
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h b/deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h
deleted file mode 100644 (file)
index 3105dde..0000000
+++ /dev/null
@@ -1,93 +0,0 @@
-#include "../include/local_assemble_plastic_project.h"
-#include <deal.II/base/quadrature_lib.h>
-
-template <int dim>
-LocalAssemblePlasticProject<dim>::LocalAssemblePlasticProject() : 
-fe(0, "Local Assemble Fe Pointer"),
-  fe_v(0, "Local Assemble FeValues Pointer")
-{
-}
-
-
-template <int dim>
-void LocalAssemblePlasticProject<dim>::reinit (FiniteElement<dim> &myfe,
-                                              Table<2, SymmetricTensor<2,dim> > &my_plastic_strain) 
-{ 
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  fe = &myfe;
-  QGauss<dim> quadrature(2*fe->degree + 1);
-  QGauss<dim-1> face_quadrature(2*fe->degree + 1);
-  UpdateFlags flags (update_gradients |
-                    update_q_points  |
-                    update_JxW_values);
-  fe_v = new FEValues<dim>(*fe, quadrature, flags);
-    
-  this->flags = assemble_rhs_cell;
-    
-  plastic_strain = my_plastic_strain;
-}
-
-template <int dim>
-LocalAssemblePlasticProject<dim>::~LocalAssemblePlasticProject()
-{
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  fe = 0;
-}
-
-template <int dim>
-void LocalAssemblePlasticProject<dim>::parameters(ParameterHandler &prm)
-{
-  prm.enter_subsection("Elastic Moduli");
-  C.parse_parameters(prm);
-  prm.leave_subsection();
-}
-
-
-template <int dim>
-void LocalAssemblePlasticProject<dim>::assemble_rhs_term
-(const typename MGDoFHandler<dim>::active_cell_iterator& cell, 
- Vector<double> &cell_rhs)
-{
-  cell_rhs = 0;
-  Assert(fe, ExcNotInitialized());
-
-  fe_v->reinit(cell);
-  
-  unsigned int size = fe->n_components();
-  unsigned int n_q_points = fe_v->n_quadrature_points;
-  
-  std::vector<Vector<double> > load_vector (n_q_points, Vector<double>(size) );
-  std::vector<Point<dim> > points = fe_v->get_quadrature_points();  
-
-  for (unsigned int qp=0; qp<n_q_points; ++qp) {
-    
-    SymmetricTensor<4,dim> C_qp = C(points[qp]);
-
-    for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
-      unsigned int comp_i = fe->system_to_component_index(i).first;
-      
-      for(unsigned int b=0; b<dim; ++b) {
-       for(unsigned int m=0; m<dim; ++m) {
-         for(unsigned int n=0; n<dim; ++n) {
-          
-           cell_rhs(i) += ( C_qp[comp_i][b][m][n] *
-                            plastic_strain(cell->index(), qp)[m][n] *
-                            fe_v->shape_grad(i, qp)[b] *
-                            fe_v->JxW(qp) );
-           
-         } //n
-       } //m 
-      } //b
-    } //i
-  } //qp
-}
-
-
diff --git a/deal.II/examples/step-41/include/local_assemble_scalar_project.h b/deal.II/examples/step-41/include/local_assemble_scalar_project.h
deleted file mode 100644 (file)
index 6c05742..0000000
+++ /dev/null
@@ -1,51 +0,0 @@
-#ifndef LOCAL_ASSEMBLE_SCALAR_PROJECT
-#define LOCAL_ASSEMBLE_SCALAR_PROJECT
-
-#include "local_assemble.h"
-#include <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/base/parameter_handler.h>
-
-
-template <int dim>
-class LocalAssembleScalarProject : public LocalAssembleBase<dim>
-{
- public:
-  
-  LocalAssembleScalarProject(); 
-  
-  ~LocalAssembleScalarProject(); 
-  
-  virtual void assemble_rhs_term
-    (const typename MGDoFHandler<dim>::active_cell_iterator&,
-     Vector<double> &);
-
-
-  void reinit(FiniteElement<dim>&,
-             Table<2, double > &,
-             Table<2, double > &);
-                       
- private:
-  /** A pointer to the finite element.*/
-  SmartPointer<FiniteElement<dim> > fe;
-
-  /** A pointer to the plastic strain */
-  Table<dim,double >   hard_table;
-  Table<dim,double >   iter_table;
-
-  /** A pointer to fe_values objects. */
-  SmartPointer<FEValues<dim> > fe_v;
-
-  
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h b/deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h
deleted file mode 100644 (file)
index 3821b90..0000000
+++ /dev/null
@@ -1,84 +0,0 @@
-#include "../include/local_assemble_scalar_project.h"
-#include <deal.II/base/quadrature_lib.h>
-
-template <int dim>
-LocalAssembleScalarProject<dim>::LocalAssembleScalarProject() : 
-fe(0, "Local Assemble Fe Pointer"),
-  fe_v(0, "Local Assemble FeValues Pointer")
-{
-}
-
-
-template <int dim>
-void LocalAssembleScalarProject<dim>::reinit (FiniteElement<dim> &myfe,
-                                             Table<2, double> &my_hard,
-                                             Table<2, double> &my_iter) 
-{ 
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  fe = &myfe;
-  QGauss<dim> quadrature(2*fe->degree + 1);
-  QGauss<dim-1> face_quadrature(2*fe->degree + 1);
-  UpdateFlags flags (update_values |
-                    update_q_points  |
-                    update_JxW_values);
-  fe_v = new FEValues<dim>(*fe, quadrature, flags);
-    
-  this->flags = assemble_rhs_cell;
-    
-  hard_table = my_hard;
-  iter_table = my_iter;
-}
-
-template <int dim>
-LocalAssembleScalarProject<dim>::~LocalAssembleScalarProject()
-{
-  if(fe_v) {
-    FEValues<dim> * p = fe_v;
-    fe_v = 0;
-    delete p;
-  }
-  fe = 0;
-}
-
-
-template <int dim>
-void LocalAssembleScalarProject<dim>::assemble_rhs_term
-(const typename MGDoFHandler<dim>::active_cell_iterator& cell, 
- Vector<double> &cell_rhs)
-{
-  cell_rhs = 0;
-  Assert(fe, ExcNotInitialized());
-
-  fe_v->reinit(cell);
-  
-  unsigned int n_q_points = fe_v->n_quadrature_points; 
-
-  for (unsigned int qp=0; qp<n_q_points; ++qp) {
-
-    for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
-      unsigned int comp_i = fe->system_to_component_index(i).first;
-      
-      if(comp_i == 0) {
-       
-       cell_rhs(i) += ( hard_table(cell->index(), qp) *
-                        fe_v->shape_value(i, qp) *
-                        fe_v->JxW(qp) );
-      }
-
-      if(comp_i == 1) {
-       
-       cell_rhs(i) += ( iter_table(cell->index(), qp) *
-                        fe_v->shape_value(i, qp) *
-                        fe_v->JxW(qp) );
-      }
-
-
-    } //i
-  } //qp
-}
-
-
diff --git a/deal.II/examples/step-41/include/my_tools.h b/deal.II/examples/step-41/include/my_tools.h
deleted file mode 100644 (file)
index dd80259..0000000
+++ /dev/null
@@ -1,989 +0,0 @@
-//---------------------------------------------------------------------------
-//    $Id: fe.h,v 1.124 2005/09/17 09:19:19 guido Exp $
-//    Version: $Name:  $
-//
-//    Copyright (C) 2005 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//---------------------------------------------------------------------------
-#ifndef __deal2__my_tools_h
-#define __deal2__my_tools_h
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/multigrid/mg_base.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/base/quadrature_lib.h>
-
-#include <deal.II/dofs/dof_constraints.h>
-#include "assemble_flags.h"
-#include "local_assemble.h"
-#include "vector_space.h"
-
-/**
- * Some assembly routines. These are an extension of the already
- * existing deal.II assembly routines that make use of the support
- * class LocalAssembleBase. These routines call functions of the local
- * assemblers passing local matrices or vectors and the cell. It is
- * responsability of the local assemblers to fill in the informations
- * correctly. Unlike step-12 of the deal.II library, the cell iterator
- * is passed, not the fe_values. This means that the local assembler
- * should build its own fe_values object. It is its responsability to
- * select all the outer details on how to actually compute these
- * matrices...
- * 
- *
- * @author Luca Heltai,
- * 2005, 2008
- */
-
-namespace MyTools {
-    using namespace std;
-
-    /** DG-Like assembly routine for matrices. Hanging node
-     * constraints are not taken into account, and one should do this
-     * outside. However face terms are correctly handled, and this
-     * method is capable of assemblying fully DG objects. */
-    template <int dim, typename DH, typename MATRIX>
-    void assemble(DH &dof_handler, MATRIX& system_matrix, 
-                 LocalAssembleBase<dim, DH> &local, 
-                 const unsigned int this_mpi_process = 0) 
-    {
-       const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
-       std::vector<unsigned int> dofs (dofs_per_cell);
-       std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
-    
-       // Now we create the cell matrices and vectors. Here we need two
-       // cell matrices, both for face terms that include test functions
-       // <code>vi</code> (internal shape functions, i.e. shape functions
-       // of the current cell). 
-       // 
-       // i stands for internal, e stands for external
-       FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
-       FullMatrix<double> ui_ve_matrix (dofs_per_cell, dofs_per_cell);
-
-       // Furthermore we need some cell iterators.
-       typename DH::active_cell_iterator
-           cell = dof_handler.begin_active(),
-           endc = dof_handler.end();
-
-       // Now we start the loop over all
-       // active cells.
-       for (;cell!=endc; ++cell) 
-           if( cell->subdomain_id() == this_mpi_process )
-           {
-               // On each cell we need to reset the <code>ui_vi_matrix</code>
-               // and <code>cell_vector</code> to zero, before assembling the
-               // cell terms.
-               ui_vi_matrix = 0;
-      
-               // and call the function that assembles the cell terms. The
-               // first argument is the <code>FEValues</code> that was
-               // previously reinit'ed on the current cell.
-               if( local.flags & assemble_cell ) 
-                   local.assemble_cell_term(cell, ui_vi_matrix);
-         
-               // As in previous examples the vector `dofs' includes the
-               // dof_indices.
-               cell->get_dof_indices (dofs);
-         
-               // This is the start of the nested loop over all faces. 
-               if( local.flags & (assemble_boundary|assemble_face) ) {
-                   for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-                   {
-                       // First we set the face iterator
-                       typename DH::face_iterator face=cell->face(face_no);
-         
-                       // and clear the <code>ui_ve_matrix</code> on each face.
-                       ui_ve_matrix = 0;
-                 
-                       // Now we distinguish the four different cases in the
-                       // ordering mentioned above. We start with faces belonging
-                       // to the boundary of the domain.
-                       if (face->at_boundary())
-                       {
-                           // and assemble the corresponding face terms.
-                           if(local.flags & assemble_boundary) 
-                               local.assemble_boundary_term(cell, face_no,
-                                                            ui_vi_matrix);
-                       }
-                       else if(local.flags & assemble_face) 
-                       {
-                           // Now we are not on the boundary of the domain,
-                           // therefore there must exist a neighboring cell.
-                           typename DH::cell_iterator neighbor=
-                               cell->neighbor(face_no);
-
-                           // We proceed with the second and most complicated case:
-                           // the neighboring cell is more refined than the current
-                           // cell. As in deal.II neighboring cells are restricted
-                           // to have a level difference of not more than one, the
-                           // neighboring cell is known to be at most ONCE more
-                           // refined than the current cell. Furthermore also the
-                           // face is more refined, i.e. it has children. Here we
-                           // note that the following part of code will not work
-                           // for <code>dim==1</code>.
-                           if (face->has_children())
-                           {
-                               // First we store which number the current cell has
-                               // in the list of neighbors of the neighboring
-                               // cell. Hence,
-                               // neighbor-@>neighbor(neighbor_face_no) equals the
-                               // current cell <code>cell</code>.
-                               const unsigned int neighbor_face_no=
-                                   cell->neighbor_of_neighbor(face_no);
-                 
-                 
-                               // We loop over subfaces
-                               for (unsigned int subface_no=0;
-                                    subface_no<face->n_children(); ++subface_no)
-                               {
-                                   // and set the cell iterator
-                                   // <code>neighbor_child</code> to the cell
-                                   // placed `behind' the current subface.
-                                   typename DH::active_cell_iterator
-                                       neighbor_child
-                                       = cell->neighbor_child_on_subface (face_no, subface_no);
-                     
-                                   // As these are quite complicated indirections
-                                   // which one does not usually get right at first
-                                   // attempt we check for the internal
-                                   // consistency.
-                                   Assert (neighbor_child->face(neighbor_face_no) == face->child(subface_no),
-                                           ExcInternalError());
-                                   Assert (!neighbor_child->has_children(), ExcInternalError());
-
-                                   // We need to reset the
-                                   // <code>ui_ve_matrix</code> on each subface
-                                   // because on each subface the <code>un</code>
-                                   // belong to different neighboring cells.
-                                   ui_ve_matrix = 0;
-                     
-                                   // As already mentioned above for the current
-                                   // case (case 2) we employ the
-                                   // <code>FESubfaceValues</code> of the current
-                                   // cell (here reinited for the current cell,
-                                   // face and subface) and we employ the
-                                   // FEFaceValues of the neighboring child cell.
-                                   local.assemble_face_term(cell, face_no, subface_no,
-                                                            neighbor_child, neighbor_face_no,
-                                                            ui_vi_matrix,
-                                                            ui_ve_matrix);
-                     
-                                   // Then we get the dof indices of the
-                                   // neighbor_child cell
-                                   neighbor_child->get_dof_indices (dofs_neighbor);
-                                                               
-                                   // and distribute <code>ui_ve_matrix</code> to
-                                   // the system_matrix
-                                   for (unsigned int i=0; i<dofs_per_cell; ++i)
-                                       for (unsigned int k=0; k<dofs_per_cell; ++k)
-                                           system_matrix.add(dofs_neighbor[i], dofs[k],
-                                                             ui_ve_matrix(i,k));
-                               }
-                               // End of <code>if (face-@>has_children())</code>
-                           }
-                           else
-                           {
-                               // We proceed with case 3, i.e. neighboring cell is
-                               // of the same refinement level as the current cell.
-                               if (neighbor->level() == cell->level()) 
-                               {
-                                   // Like before we store which number the current
-                                   // cell has in the list of neighbors of the
-                                   // neighboring cell.
-                                   const unsigned int neighbor_face_no=cell->neighbor_of_neighbor(face_no);
-
-                                   local.assemble_face_term(cell, face_no,
-                                                            neighbor, neighbor_face_no,
-                                                            ui_vi_matrix,
-                                                            ui_ve_matrix);
-                                   // End of <code>if (neighbor-@>level() ==
-                                   // cell-@>level())</code>
-                               }
-                               else
-                               {
-                                   // Finally we consider case 4. When the
-                                   // neighboring cell is not finer and not of the
-                                   // same refinement level as the current cell it
-                                   // must be coarser.
-                                   Assert(neighbor->level() < cell->level(), ExcInternalError());
-
-                                   // Find out the how many'th face_no and
-                                   // subface_no the current face is w.r.t. the
-                                   // neighboring cell.
-                                   const std::pair<unsigned int, unsigned int> faceno_subfaceno=
-                                       cell->neighbor_of_coarser_neighbor(face_no);
-                                   const unsigned int neighbor_face_no=faceno_subfaceno.first,
-                                       neighbor_subface_no=faceno_subfaceno.second;
-
-                                   Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
-                                                                                neighbor_subface_no)
-                                           == cell,
-                                           ExcInternalError());
-
-                                   local.assemble_face_term(cell, face_no,
-                                                            neighbor, 
-                                                            neighbor_face_no,
-                                                            neighbor_subface_no,
-                                                            ui_vi_matrix,
-                                                            ui_ve_matrix);
-                               }
-
-                               // Now we get the dof indices of the
-                               // <code>neighbor</code> cell,
-                               neighbor->get_dof_indices (dofs_neighbor);
-                                                               
-                               // and distribute the
-                               // <code>ui_ve_matrix</code>.
-                               for (unsigned int i=0; i<dofs_per_cell; ++i)
-                                   for (unsigned int k=0; k<dofs_per_cell; ++k)
-                                       system_matrix.add(dofs_neighbor[i], dofs[k],
-                                                         ui_ve_matrix(i,k));
-                           }
-                           // End of <code>face not at boundary</code>:
-                       }
-                       // End of loop over all faces:
-                   }
-               }
-           
-               // Finally we distribute the
-               // <code>ui_vi_matrix</code>
-               for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
-           }
-    }
-
-
-    /** MG assembly routine for matrices. Hanging node constraints are
-     * taken into account. However face terms are NOT correctly
-     * handled, that is, this method is currently not capable of
-     * handling fluxex. 
-     
-     Different levels can be assembled. If the level is -1, then the
-     active level is used.
-     */
-    template <int dim, typename MATRIX>
-    void assemble(MGDoFHandler<dim> &mg_dof_handler, ConstraintMatrix &hang,
-                 int level, MATRIX& m, LocalAssembleBase<dim, MGDoFHandler<dim> > &local, 
-                 const unsigned int this_mpi_process = 0) 
-    {
-       AssertThrow( !(local.flags & assemble_face), ExcNotImplemented());
-       // Cell iterators.
-       typename MGDoFHandler<dim>::active_cell_iterator cell, endc;
-       if(level < 0) { 
-           cell = mg_dof_handler.begin_active();
-           endc = mg_dof_handler.end();
-       } else {
-           cell = mg_dof_handler.begin_active(level);  
-           endc = mg_dof_handler.end(level);
-       }
-
-       unsigned int dofs_per_cell = 
-           mg_dof_handler.get_fe().dofs_per_cell;
-
-       // dofs and neighbor dofs.
-       vector<unsigned int> dofs(dofs_per_cell);
-       vector<unsigned int> neighbor_dofs(dofs_per_cell);
-
-       FullMatrix<double> cell_m(dofs_per_cell, dofs_per_cell);
-       FullMatrix<double> cell_neighbor_m(dofs_per_cell, dofs_per_cell);
-
-       for (; cell!=endc; ++cell)
-           if( cell->subdomain_id() == this_mpi_process )
-           {
-               // Reset local matrices.
-               cell_m = 0;
-               cell_neighbor_m = 0;
-         
-               // Update indices
-               if(level<0)
-                   cell->get_dof_indices(dofs);
-               else
-                   cell->get_mg_dof_indices (dofs);
-
-               // Assemble cell term if needed
-               if( local.flags & assemble_cell ) {
-
-                 //
-                 local.assemble_cell_term(cell, cell_m);
-               }
-               for(unsigned int n=0; n<dofs.size(); ++n)
-               hang.distribute_local_to_global(cell_m, dofs, m);
-#if (DEAL_II_DIM > 1)
-               // Assemble boundary or face terms if needed
-               if( local.flags & (assemble_boundary|assemble_face) ) {
-           
-                   for (unsigned int face_no=0; 
-                        face_no<GeometryInfo<dim>::faces_per_cell; 
-                        ++face_no)
-                   {
-                       typename MGDoFHandler<dim>::face_iterator face=cell->face(face_no);
-               
-                       if (face->at_boundary())
-                       {
-                           if(local.flags & assemble_boundary)   {
-                               cell_m = 0;
-                               local.assemble_boundary_term(cell, face_no, cell_m);
-                               hang.distribute_local_to_global(cell_m, dofs, m);
-                           }
-                       } // End of Boundary terms
-                       else if(local.flags & assemble_face) {
-                           typename MGDoFHandler<dim>::cell_iterator neighbor=
-                               cell->neighbor(face_no);;
-                 
-                           if (face->has_children()) // The neighbor is more refined
-                           {
-                               const unsigned int n_face_no=
-                                   cell->neighbor_of_neighbor(face_no);
-                     
-                     
-                               for (unsigned int subface_no=0;
-                                    subface_no<face->n_children(); ++subface_no)
-                               {
-                                   typename MGDoFHandler<dim>::active_cell_iterator
-                                       neighbor_child
-                                       = cell->neighbor_child_on_subface (face_no, subface_no);
-                         
-                                   Assert (neighbor_child->face(n_face_no) == face->child(subface_no),
-                                           ExcInternalError());
-                                   Assert (!neighbor_child->has_children(), ExcInternalError());
-                         
-                                   cell_m = 0;
-                                   cell_neighbor_m = 0;
-                                   local.assemble_face_term(cell, face_no, subface_no,
-                                                            neighbor_child, n_face_no,
-                                                            cell_m, cell_neighbor_m);
-                                   if(level<0)
-                                       neighbor_child->get_dof_indices (neighbor_dofs);
-                                   else 
-                                       neighbor_child->get_mg_dof_indices (neighbor_dofs);
-                         
-                                   hang.distribute_local_to_global(cell_m, dofs, m);
-                                   // TBA Distribute neighbors!!!
-                                   // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
-                                   //
-
-
-                               }
-                           } // face has children
-                           else
-                           {
-                               if (neighbor->level() == cell->level())
-                               {
-                                   const unsigned int n_face_no=cell->neighbor_of_neighbor(face_no);
-                         
-                                   cell_neighbor_m = 0;
-                                   cell_m = 0;
-                                   local.assemble_face_term(cell, face_no, 
-                                                            neighbor, n_face_no, 
-                                                            cell_m,
-                                                            cell_neighbor_m);
-                         
-                                   if(level<0)
-                                       neighbor->get_dof_indices (neighbor_dofs);
-                                   else 
-                                       neighbor->get_mg_dof_indices (neighbor_dofs);
-                         
-                                   hang.distribute_local_to_global(cell_m, dofs, m);
-                                   // TBA: distribute neighbor dofs....
-                                   // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
-                                   // 
-                         
-                               } // Same level on the two neighbors
-                               else
-                               {
-                                   Assert(neighbor->level() < cell->level(), ExcInternalError());
-                         
-                                   const std::pair<unsigned int, unsigned int> faceno_subfaceno=
-                                       cell->neighbor_of_coarser_neighbor(face_no);
-                                   const unsigned int neighbor_face_no=faceno_subfaceno.first,
-                                       neighbor_subface_no=faceno_subfaceno.second;
-                         
-                                   Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
-                                                                                neighbor_subface_no)
-                                           == cell,
-                                           ExcInternalError());
-                         
-                                   cell_neighbor_m = 0;
-                                   cell_m = 0;
-                                   local.assemble_face_term(cell, face_no,
-                                                            neighbor, neighbor_face_no, 
-                                                            neighbor_subface_no,
-                                                            cell_m, cell_neighbor_m);
-                         
-                                   if(level<0)
-                                       neighbor->get_dof_indices (neighbor_dofs);
-                                   else 
-                                       neighbor->get_mg_dof_indices (neighbor_dofs);
-                         
-                                   hang.distribute_local_to_global(cell_m, dofs, m);
-                                   // TBA Distribute neighbor dofs
-                                   // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
-                                   //
-
-                               } // The neighbor is coarser
-                           } // else : if face has children
-                       } // skipped if no flags for face terms
-                   } // faces loop
-               } // Skipped if no flags on either faces or boundary 
-#else
-    //1D face terms
-#endif 
-           } // cells loop
-    }
-  
-    /** MG assembly routine for matrices and vectors. Hanging node
-     * constraints are taken into account. However face terms are NOT
-     * correctly handled, that is, this method is currently not
-     * capable of handling fluxex.
-     
-     Different levels can be assembled. If the level is -1, then the
-     active level is used.
-     */
-    template <int dim, typename MATRIX, typename VEC>
-    void assemble(MGDoFHandler<dim> &mg_dof_handler, ConstraintMatrix &hang,
-                 int level, MATRIX& m, VEC& rhs, LocalAssembleBase<dim, MGDoFHandler<dim> > &local, 
-                 const unsigned int this_mpi_process = 0) 
-    {
-       // Cell iterators.
-       typename MGDoFHandler<dim>::active_cell_iterator cell, endc;
-       if(level < 0) { 
-           cell = mg_dof_handler.begin_active();
-           endc = mg_dof_handler.end();
-       } else {
-           cell = mg_dof_handler.begin_active(level);  
-           endc = mg_dof_handler.end(level);
-       }
-
-       unsigned int dofs_per_cell = 
-           mg_dof_handler.get_fe().dofs_per_cell;
-
-       // dofs and neighbor dofs.
-       vector<unsigned int> dofs(dofs_per_cell);
-       vector<unsigned int> neighbor_dofs(dofs_per_cell);
-
-       FullMatrix<double> cell_m(dofs_per_cell, dofs_per_cell);
-       FullMatrix<double> cell_neighbor_m(dofs_per_cell, dofs_per_cell);
-       Vector<double>     local_rhs(dofs_per_cell);
-
-       for (; cell!=endc; ++cell)
-           if( cell->subdomain_id() == this_mpi_process )
-           {
-               // Reset local matrices.
-               cell_m = 0;
-               cell_neighbor_m = 0;
-         
-               // Update indices
-               if(level<0)
-                   cell->get_dof_indices(dofs);
-               else
-                   cell->get_mg_dof_indices (dofs);
-
-               // Assemble cell term if needed
-               if( local.flags & assemble_cell ) 
-                   local.assemble_cell_term(cell, cell_m);
-               if( local.flags & assemble_rhs_cell)
-                   local.assemble_rhs_term(cell, local_rhs);
-               
-               hang.distribute_local_to_global(cell_m, dofs, m);
-               hang.distribute_local_to_global(local_rhs, dofs, rhs);
-
-#if (DEAL_II_DIM > 1)
-               // Assemble boundary or face terms if needed
-               if( local.flags & (assemble_boundary|assemble_face) ) {
-           
-                   for (unsigned int face_no=0; 
-                        face_no<GeometryInfo<dim>::faces_per_cell; 
-                        ++face_no)
-                   {
-                       typename MGDoFHandler<dim>::face_iterator face=cell->face(face_no);
-               
-                       if (face->at_boundary())
-                       {
-                           if(local.flags & assemble_boundary)   {
-                               cell_m = 0;
-                               local.assemble_boundary_term(cell, face_no, cell_m);
-                               hang.distribute_local_to_global(cell_m, dofs, m);
-                           }
-                           if(local.flags & assemble_rhs_boundary) {
-                               local_rhs =0;
-                               local.assemble_rhs_boundary_term(cell, face_no, local_rhs);
-                               hang.distribute_local_to_global(local_rhs, dofs, rhs);
-                           }
-                       } // End of Boundary terms
-                       else if(local.flags & assemble_face) {
-                           typename MGDoFHandler<dim>::cell_iterator neighbor=
-                               cell->neighbor(face_no);;
-                 
-                           if (face->has_children()) // The neighbor is more refined
-                           {
-                               const unsigned int n_face_no=
-                                   cell->neighbor_of_neighbor(face_no);
-                     
-                     
-                               for (unsigned int subface_no=0;
-                                    subface_no<face->n_children(); ++subface_no)
-                               {
-                                   typename MGDoFHandler<dim>::active_cell_iterator
-                                       neighbor_child
-                                       = cell->neighbor_child_on_subface (face_no, subface_no);
-                         
-                                   Assert (neighbor_child->face(n_face_no) == face->child(subface_no),
-                                           ExcInternalError());
-                                   Assert (!neighbor_child->has_children(), ExcInternalError());
-                         
-                                   cell_m = 0;
-                                   cell_neighbor_m = 0;
-                                   local.assemble_face_term(cell, face_no, subface_no,
-                                                            neighbor_child, n_face_no,
-                                                            cell_m, cell_neighbor_m);
-                                   if(level<0)
-                                       neighbor_child->get_dof_indices (neighbor_dofs);
-                                   else 
-                                       neighbor_child->get_mg_dof_indices (neighbor_dofs);
-                         
-                                   hang.distribute_local_to_global(cell_m, dofs, m);
-                                   // TBA Distribute neighbors!!!
-                                   // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
-                                   //
-
-
-                               }
-                           } // face has children
-                           else
-                           {
-                               if (neighbor->level() == cell->level())
-                               {
-                                   const unsigned int n_face_no=cell->neighbor_of_neighbor(face_no);
-                         
-                                   cell_neighbor_m = 0;
-                                   cell_m = 0;
-                                   local.assemble_face_term(cell, face_no, 
-                                                            neighbor, n_face_no, 
-                                                            cell_m,
-                                                            cell_neighbor_m);
-                         
-                                   if(level<0)
-                                       neighbor->get_dof_indices (neighbor_dofs);
-                                   else 
-                                       neighbor->get_mg_dof_indices (neighbor_dofs);
-                         
-                                   hang.distribute_local_to_global(cell_m, dofs, m);
-                                   // TBA: distribute neighbor dofs....
-                                   // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
-                                   // 
-                         
-                               } // Same level on the two neighbors
-                               else
-                               {
-                                   Assert(neighbor->level() < cell->level(), ExcInternalError());
-                         
-                                   const std::pair<unsigned int, unsigned int> faceno_subfaceno=
-                                       cell->neighbor_of_coarser_neighbor(face_no);
-                                   const unsigned int neighbor_face_no=faceno_subfaceno.first,
-                                       neighbor_subface_no=faceno_subfaceno.second;
-                         
-                                   Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
-                                                                                neighbor_subface_no)
-                                           == cell,
-                                           ExcInternalError());
-                         
-                                   cell_neighbor_m = 0;
-                                   cell_m = 0;
-                                   local.assemble_face_term(cell, face_no,
-                                                            neighbor, neighbor_face_no, 
-                                                            neighbor_subface_no,
-                                                            cell_m, cell_neighbor_m);
-                         
-                                   if(level<0)
-                                       neighbor->get_dof_indices (neighbor_dofs);
-                                   else 
-                                       neighbor->get_mg_dof_indices (neighbor_dofs);
-                         
-                                   hang.distribute_local_to_global(cell_m, dofs, m);
-                                   // TBA Distribute neighbor dofs
-                                   // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
-                                   //
-
-                               } // The neighbor is coarser
-                           } // else : if face has children
-                       } // skipped if no flags for face terms
-                   } // faces loop
-
-               } // Skipped if no flags on either faces or boundary 
-
-#else
-    //1D face terms
-#endif 
-           } // cells loop
-    }
-
-
-      /** MG assembly routine for rhs vectors. Hanging node
-     * constraints are taken into account. However face terms are NOT
-     * correctly handled, that is, this method is currently not
-     * capable of handling fluxex.
-     
-     Different levels can be assembled. If the level is -1, then the
-     active level is used.
-     */
-    template <int dim, typename VECTOR>
-    void assemble_rhs(MGDoFHandler<dim> &mg_dof_handler, ConstraintMatrix &hang,
-                     int level, VECTOR& rhs, LocalAssembleBase<dim, MGDoFHandler<dim> > &local,
-                     const unsigned int this_mpi_process = 0) 
-    {
-       // Cell iterators.
-       typename MGDoFHandler<dim>::active_cell_iterator cell, endc;
-       if(level < 0) { 
-           cell = mg_dof_handler.begin_active();
-           endc = mg_dof_handler.end();
-       } else {
-           cell = mg_dof_handler.begin_active(level);  
-           endc = mg_dof_handler.end(level);
-       }
-      
-       unsigned int dofs_per_cell = 
-           mg_dof_handler.get_fe().dofs_per_cell;
-      
-       // dofs and neighbor dofs.
-       vector<unsigned int> dofs(dofs_per_cell);
-      
-       Vector<double> cell_rhs(dofs_per_cell);
-      
-       for (; cell!=endc; ++cell)
-           if( cell->subdomain_id() == this_mpi_process )
-           {
-               // Update indices
-               if(level<0)
-                   cell->get_dof_indices(dofs);
-               else
-                   cell->get_mg_dof_indices (dofs);
-         
-               // Assemble cell term if needed
-               if( local.flags & assemble_rhs_cell ) {
-                   // Reset the cell term
-                   cell_rhs = 0;
-                   local.assemble_rhs_term(cell, cell_rhs);
-                   // Distribute the local rhs to the global matrix.
-                   hang.distribute_local_to_global(cell_rhs, dofs, rhs);
-               }
-
-               // Assemble boundary or face terms if needed
-               if( local.flags & (assemble_rhs_boundary/* |assemble_rhs_face */) )
-                   for (unsigned int face_no=0; 
-                        face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-                   {
-                       typename MGDoFHandler<dim>::face_iterator face=cell->face(face_no);
-                       if (face->at_boundary())
-                       { 
-                           cell_rhs = 0;
-                           local.assemble_rhs_boundary_term(cell, face_no, cell_rhs);
-                           // Distribute the local rhs to the global matrix.
-                           hang.distribute_local_to_global(cell_rhs, dofs, rhs);
-                       }
-                   } // face loop        
-           } // Cell loop
-    }
-    
-    /** DG-Like assembly routine for matrices and vectors. Hanging
-     * node constraints are not taken into account, and one should do
-     * this outside. However face terms are correctly handled, and
-     * this method is capable of assemblying fully DG objects. */
-    template <int dim, typename DH, typename MATRIX, typename VEC>
-      void assemble(DH &dof_handler, MATRIX& system_matrix, VEC& system_rhs,
-                 LocalAssembleBase<dim, DH> &local, 
-                 const unsigned int this_mpi_process = 0) 
-    {
-       const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
-       std::vector<unsigned int> dofs (dofs_per_cell);
-       std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
-    
-       // Now we create the cell matrices and vectors. Here we need two
-       // cell matrices, both for face terms that include test functions
-       // <code>vi</code> (internal shape functions, i.e. shape functions
-       // of the current cell). 
-       // 
-       // i stands for internal, e stands for external
-       FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
-       FullMatrix<double> ui_ve_matrix (dofs_per_cell, dofs_per_cell);
-       Vector<double>     local_rhs(dofs_per_cell);
-
-       // Furthermore we need some cell iterators.
-       typename DH::active_cell_iterator
-           cell = dof_handler.begin_active(),
-           endc = dof_handler.end();
-
-       // Now we start the loop over all
-       // active cells.
-       for (;cell!=endc; ++cell) 
-           if( cell->subdomain_id() == this_mpi_process )
-           {
-               // On each cell we need to reset the <code>ui_vi_matrix</code>
-               // and <code>cell_vector</code> to zero, before assembling the
-               // cell terms.
-               ui_vi_matrix = 0;
-               local_rhs = 0;
-      
-               // and call the function that assembles the cell terms. The
-               // first argument is the <code>FEValues</code> that was
-               // previously reinit'ed on the current cell.
-               if( local.flags & assemble_cell ) 
-                   local.assemble_cell_term(cell, ui_vi_matrix);
-               
-               if( local.flags & assemble_rhs_cell )
-                   local.assemble_rhs_term(cell, local_rhs);
-         
-               // As in previous examples the vector `dofs' includes the
-               // dof_indices.
-               cell->get_dof_indices (dofs);
-         
-               // This is the start of the nested loop over all faces. 
-               if( local.flags & (assemble_boundary|assemble_face) ) {
-                   for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-                   {
-                       // First we set the face iterator
-                       typename DH::face_iterator face=cell->face(face_no);
-         
-                       // and clear the <code>ui_ve_matrix</code> on each face.
-                       ui_ve_matrix = 0;
-                 
-                       // Now we distinguish the four different cases in the
-                       // ordering mentioned above. We start with faces belonging
-                       // to the boundary of the domain.
-                       if (face->at_boundary())
-                       {
-                           // and assemble the corresponding face terms.
-                           if(local.flags & assemble_boundary) 
-                               local.assemble_boundary_term(cell, face_no,
-                                                            ui_vi_matrix);
-                           if(local.flags & assemble_rhs_boundary)
-                               local.assemble_rhs_boundary_term(cell, face_no,
-                                                                local_rhs);
-                       }
-                       else if(local.flags & assemble_face) 
-                       {
-                           // Now we are not on the boundary of the domain,
-                           // therefore there must exist a neighboring cell.
-                           typename DH::cell_iterator neighbor=
-                               cell->neighbor(face_no);
-
-                           // We proceed with the second and most complicated case:
-                           // the neighboring cell is more refined than the current
-                           // cell. As in deal.II neighboring cells are restricted
-                           // to have a level difference of not more than one, the
-                           // neighboring cell is known to be at most ONCE more
-                           // refined than the current cell. Furthermore also the
-                           // face is more refined, i.e. it has children. Here we
-                           // note that the following part of code will not work
-                           // for <code>dim==1</code>.
-                           if (face->has_children())
-                           {
-                               // First we store which number the current cell has
-                               // in the list of neighbors of the neighboring
-                               // cell. Hence,
-                               // neighbor-@>neighbor(neighbor_face_no) equals the
-                               // current cell <code>cell</code>.
-                               const unsigned int neighbor_face_no=
-                                   cell->neighbor_of_neighbor(face_no);
-                 
-                 
-                               // We loop over subfaces
-                               for (unsigned int subface_no=0;
-                                    subface_no<face->n_children(); ++subface_no)
-                               {
-                                   // and set the cell iterator
-                                   // <code>neighbor_child</code> to the cell
-                                   // placed `behind' the current subface.
-                                   typename DH::active_cell_iterator
-                                       neighbor_child
-                                       = cell->neighbor_child_on_subface (face_no, subface_no);
-                     
-                                   // As these are quite complicated indirections
-                                   // which one does not usually get right at first
-                                   // attempt we check for the internal
-                                   // consistency.
-                                   Assert (neighbor_child->face(neighbor_face_no) == face->child(subface_no),
-                                           ExcInternalError());
-                                   Assert (!neighbor_child->has_children(), ExcInternalError());
-
-                                   // We need to reset the
-                                   // <code>ui_ve_matrix</code> on each subface
-                                   // because on each subface the <code>un</code>
-                                   // belong to different neighboring cells.
-                                   ui_ve_matrix = 0;
-                     
-                                   // As already mentioned above for the current
-                                   // case (case 2) we employ the
-                                   // <code>FESubfaceValues</code> of the current
-                                   // cell (here reinited for the current cell,
-                                   // face and subface) and we employ the
-                                   // FEFaceValues of the neighboring child cell.
-                                   local.assemble_face_term(cell, face_no, subface_no,
-                                                            neighbor_child, neighbor_face_no,
-                                                            ui_vi_matrix,
-                                                            ui_ve_matrix);
-                     
-                                   // Then we get the dof indices of the
-                                   // neighbor_child cell
-                                   neighbor_child->get_dof_indices (dofs_neighbor);
-                                                               
-                                   // and distribute <code>ui_ve_matrix</code> to
-                                   // the system_matrix
-                                   for (unsigned int i=0; i<dofs_per_cell; ++i)
-                                       for (unsigned int k=0; k<dofs_per_cell; ++k)
-                                           system_matrix.add(dofs_neighbor[i], dofs[k],
-                                                             ui_ve_matrix(i,k));
-                               }
-                               // End of <code>if (face-@>has_children())</code>
-                           }
-                           else
-                           {
-                               // We proceed with case 3, i.e. neighboring cell is
-                               // of the same refinement level as the current cell.
-                               if (neighbor->level() == cell->level()) 
-                               {
-                                   // Like before we store which number the current
-                                   // cell has in the list of neighbors of the
-                                   // neighboring cell.
-                                   const unsigned int neighbor_face_no=cell->neighbor_of_neighbor(face_no);
-
-                                   local.assemble_face_term(cell, face_no,
-                                                            neighbor, neighbor_face_no,
-                                                            ui_vi_matrix,
-                                                            ui_ve_matrix);
-                                   // End of <code>if (neighbor-@>level() ==
-                                   // cell-@>level())</code>
-                               }
-                               else
-                               {
-                                   // Finally we consider case 4. When the
-                                   // neighboring cell is not finer and not of the
-                                   // same refinement level as the current cell it
-                                   // must be coarser.
-                                   Assert(neighbor->level() < cell->level(), ExcInternalError());
-
-                                   // Find out the how many'th face_no and
-                                   // subface_no the current face is w.r.t. the
-                                   // neighboring cell.
-                                   const std::pair<unsigned int, unsigned int> faceno_subfaceno=
-                                       cell->neighbor_of_coarser_neighbor(face_no);
-                                   const unsigned int neighbor_face_no=faceno_subfaceno.first,
-                                       neighbor_subface_no=faceno_subfaceno.second;
-
-                                   Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
-                                                                                neighbor_subface_no)
-                                           == cell,
-                                           ExcInternalError());
-
-                                   local.assemble_face_term(cell, face_no,
-                                                            neighbor, 
-                                                            neighbor_face_no,
-                                                            neighbor_subface_no,
-                                                            ui_vi_matrix,
-                                                            ui_ve_matrix);
-                               }
-
-                               // Now we get the dof indices of the
-                               // <code>neighbor</code> cell,
-                               neighbor->get_dof_indices (dofs_neighbor);
-                                                               
-                               // and distribute the
-                               // <code>ui_ve_matrix</code>.
-                               for (unsigned int i=0; i<dofs_per_cell; ++i)
-                                   for (unsigned int k=0; k<dofs_per_cell; ++k)
-                                       system_matrix.add(dofs_neighbor[i], dofs[k],
-                                                         ui_ve_matrix(i,k));
-                           }
-                           // End of <code>face not at boundary</code>:
-                       }
-                       // End of loop over all faces:
-                   }
-               }
-           
-               // Finally we distribute the
-               // <code>ui_vi_matrix</code> 
-               for (unsigned int i=0; i<dofs_per_cell; ++i) {
-                   for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
-                   system_rhs(dofs[i]) += local_rhs(i);
-               }
-           }
-    }
-    
-    /** Assembly routine for rhs. Hanging node constraints are not
-     * taken into account, and one should do this outside. */
-    template <int dim, typename DH, typename VEC>
-      void assemble_rhs(DH &dof_handler, VEC& system_rhs,
-                     LocalAssembleBase<dim, DH> &local, 
-                     const unsigned int this_mpi_process = 0) 
-    {
-       const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
-       std::vector<unsigned int> dofs (dofs_per_cell);
-       
-       Vector<double>     local_rhs(dofs_per_cell);
-
-       // Furthermore we need some cell iterators.
-       typename DH::active_cell_iterator
-           cell = dof_handler.begin_active(),
-           endc = dof_handler.end();
-
-       // Now we start the loop over all
-       // active cells.
-       for (;cell!=endc; ++cell) 
-           if( cell->subdomain_id() == this_mpi_process )
-           {
-               // On each cell we need to reset the <code>ui_vi_matrix</code>
-               // and <code>cell_vector</code> to zero, before assembling the
-               // cell terms.
-               local_rhs = 0;
-      
-               // and call the function that assembles the cell terms. The
-               // first argument is the <code>FEValues</code> that was
-               // previously reinit'ed on the current cell.
-               if( local.flags & assemble_rhs_cell )
-                   local.assemble_rhs_term(cell, local_rhs);
-         
-               // As in previous examples the vector `dofs' includes the
-               // dof_indices.
-               cell->get_dof_indices (dofs);
-         
-               // This is the start of the nested loop over all faces. 
-               if( local.flags & assemble_rhs_boundary ) {
-                   for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-                   {
-                       // First we set the face iterator
-                       typename DH::face_iterator face=cell->face(face_no);
-                 
-                       // Now we distinguish the four different cases in the
-                       // ordering mentioned above. We start with faces belonging
-                       // to the boundary of the domain.
-                       if (face->at_boundary())
-                       {
-                           if(local.flags)
-                               local.assemble_rhs_boundary_term(cell, face_no,
-                                                                local_rhs);
-                       }
-                   }
-               }
-           
-               // Finally we distribute the
-               // <code>rhs</code> 
-               for (unsigned int i=0; i<dofs_per_cell; ++i) {
-                   system_rhs(dofs[i]) += local_rhs(i);
-               }
-           }
-    }
-}
-
-#endif
diff --git a/deal.II/examples/step-41/include/output_processor.h b/deal.II/examples/step-41/include/output_processor.h
deleted file mode 100644 (file)
index 7be1dd4..0000000
+++ /dev/null
@@ -1,176 +0,0 @@
-#ifndef LH_OUTPUT_PROCESSOR_H
-#define LH_OUTPUT_PROCESSOR_H
-
-#include <fstream>
-
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_constraints.h>
-#include <deal.II/lac/block_vector.h>
-
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/filtered_iterator.h>
-
-// #include <deal.II/numerics/error_estimator.h>
-#include <deal.II/base/function.h>
-#include <deal.II/numerics/solution_transfer.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <deal.II/base/convergence_table.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/parameter_handler.h>
-
-#include "vector_space.h"
-#include <map>
-
-template<int dim, typename DH=DoFHandler<dim> >
-class FilteredDataOut : public DataOut<dim, DH>
-{
-public:
-  FilteredDataOut (const unsigned int subdomain_id)
-    :
-    subdomain_id (subdomain_id)
-  {}
-       
-  virtual typename DH::cell_iterator
-  first_cell ()
-  {
-    typename DH::active_cell_iterator
-      cell = this->dofs->begin_active();
-    while ((cell != this->dofs->end()) &&
-          (cell->subdomain_id() != subdomain_id))
-      ++cell;
-           
-    return cell;
-  }
-       
-  virtual typename DH::cell_iterator
-  next_cell (const typename DH::cell_iterator &old_cell)
-  {
-    if (old_cell != this->dofs->end())
-      {
-       const IteratorFilters::SubdomainEqualTo
-         predicate(subdomain_id);
-               
-       return
-         ++(FilteredIterator
-            <typename DH::active_cell_iterator>
-            (predicate,old_cell));
-      }
-    else
-      return old_cell;
-  }
-       
-private:
-  const unsigned int subdomain_id;
-};
-
-
-
-template <int dim, typename VECTOR=BlockVector<double> >
-class OutputProcessor : public Subscriptor
-{
-public:
-  /** The constructor takes the mpi initialization stuff. */
-  OutputProcessor (const unsigned int n_mpi_processes = 1, 
-                const unsigned int this_mpi_process = 0);
-
-  /** Initialize the given values for the paramter file. */
-  static void declare_parameters(ParameterHandler &prm);
-
-  /** Parse the given parameter handler. */
-  void parse_parameters(ParameterHandler &prm);
-  
-  /** Calculate the error of the numeric solution in variuous norms. Store
-      the result in the given table. */
-  void error_from_exact(const VectorSpace<dim> & vspace, 
-                       const VECTOR &solution, 
-                       const Function<dim> &exact,
-                       unsigned int table_no = 0,
-                       double dt=0.); 
-  
-  /** Difference between two solutions in two different vector spaces. */
-  void difference(const VectorSpace<dim> &, const VECTOR &,
-                 const VectorSpace<dim> &, const VECTOR &, 
-                 unsigned int table_no = 0, double dt=0.); 
-  
-  /** Difference between two solutions in the same vector space. */
-  void difference(const VectorSpace<dim> &, const VECTOR &, 
-                 const VECTOR &, unsigned int table_no = 0, double dt=0.); 
-  
-  /** Output solution. Using some standard format. This is now deprecated. */
-  void output_solution(const VectorSpace<dim> &, const VECTOR &, const std::string &);
-  
-  /** Prepare to output data on the given file. This will initialize
-      the data_out object and a file with a filename that is the
-      combination of the @p filename, eventually a
-      processor number and the output suffix.  */
-  void prepare_data_output(const DoFHandler<dim> &dh, const std::string &filename);
-
-  /** Add the given vector to the output file. Prior to calling this
-      method, you have to call the prepare_data_output method. The
-      string can be a comma separated list of components, or a single
-      description. In this latter case, a progressive number per
-      component is added in the end. */
-  void add_data_vector(const VECTOR &data_vector, const std::string &desc);
-  
-  /** Actually write the file. Once the data_out has been prepared,
-      vectors have been added, the data can be written to a file. This
-      is done in this class. At the end of this function call,
-      data_out and output_file are in a pristine situation, and the
-      process can be started again.*/
-  void write_data_and_clear();
-  
-  /** Dump vector. Just dump the vector to a file. This is useful to
-      read back in the future and restart from where we left off.*/
-  void dump_vector(const VECTOR &, const std::string &);
-  
-  /** By default output first table. */
-  void output_table(const unsigned int table_no=0);
-  
-private:  
-  /** Error results.*/
-  std::vector<ConvergenceTable>  tables;  
-  
-  /** Headers for tables and output. Contains the name of the solution
-      components. */
-  std::vector<std::string> headers;
-  
-  /** Headers for latex tables. Contains the name of the solution
-      components. */
-  std::vector<std::string> latex_headers;
-  
-  /** Captions for latex. */
-  std::vector<std::string> latex_captions;
-  
-  /** Names of the tables. */
-  std::vector<std::string> names;
-    
-  /** The parameters have been read. */
-  bool initialized;
-  
-  /** Write the solution. */
-  bool write_solution;
-  
-  /** Solution format. */
-  std::string solution_format;
-  
-  /** Output the partitioning of the domain. */
-  bool output_partitioning;
-  
-  /** Number of MPI processes. */
-  const unsigned int n_mpi_processes;
-  
-  /** Id of MPI process. */
-  const unsigned int this_mpi_process;
-  
-  /** Output file. */
-  std::ofstream output_file;
-  
-  /** Outputs only the data that refers to this process. */
-  FilteredDataOut<dim, DoFHandler<dim> > data_out;
-};
-
-#endif
diff --git a/deal.II/examples/step-41/include/output_processor.templates.h b/deal.II/examples/step-41/include/output_processor.templates.h
deleted file mode 100644 (file)
index daabff1..0000000
+++ /dev/null
@@ -1,135 +0,0 @@
-#include "../include/output_processor.h"
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/grid/grid_tools.h>
-
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/fe/fe.h>
-
-template <int dim, typename VECTOR>
-OutputProcessor<dim,VECTOR>::OutputProcessor (const unsigned int n_mpi, 
-                                              const unsigned int this_mpi) :
-  n_mpi_processes(n_mpi),
-  this_mpi_process(this_mpi),
-  data_out(this_mpi)
-{
-  initialized = false;
-}
-
-template <int dim, typename VECTOR>
-void OutputProcessor<dim,VECTOR>::declare_parameters (ParameterHandler &prm)
-{
-  prm.declare_entry ("Solution names", "u", Patterns::Anything(), 
-                    "This is used to output the solution to a file.");
-  
-  prm.declare_entry ("Output partitioning", "false", Patterns::Bool());
-  
-  prm.enter_subsection("Solution Out Format");
-  DataOut<dim>::declare_parameters(prm);
-  prm.set("Output format", "vtk");
-  prm.leave_subsection();
-}
-template <int dim, typename VECTOR>
-void OutputProcessor<dim,VECTOR>::parse_parameters (ParameterHandler &prm)
-{
-  output_partitioning = prm.get_bool ("Output partitioning");
-
-  headers = Utilities::split_string_list(prm.get ("Solution names"));
-
-  prm.enter_subsection("Solution Out Format");
-  data_out.parse_parameters(prm);
-  prm.leave_subsection();
-
-  initialized = true;
-}
-
-
-template <int dim, typename VECTOR>
-void OutputProcessor<dim,VECTOR>::prepare_data_output(const DoFHandler<dim> &dh,
-                                                   const std::string &filename) {
-  AssertThrow(initialized, ExcNotInitialized());
-  if(data_out.default_suffix() != "") {
-    
-      // If the output is needed and we have many processes, just output
-      // the one we need *in intermediate format*.
-      std::string fname = filename;
-      if(n_mpi_processes > 1) {
-         fname += ("_" + Utilities::int_to_string(this_mpi_process, 2) + 
-                   data_out.default_suffix(DataOutBase::deal_II_intermediate)) ;
-      } else {
-         fname += data_out.default_suffix();
-      }
-      
-      deallog << "Will write on file: " << fname.c_str() << std::endl;
-      output_file.open(fname.c_str());
-      AssertThrow(output_file, ExcIO());
-      data_out.attach_dof_handler (dh);
-      
-      if(n_mpi_processes > 1) {
-         // Output the partitioning
-         if(output_partitioning) {
-             std::vector<unsigned int> partition_int (dh.get_tria().n_active_cells());
-             GridTools::get_subdomain_association (dh.get_tria(), partition_int);
-             Vector<double> partitioning(partition_int.begin(),
-                                         partition_int.end());
-             static Vector<double> static_partitioning;
-             static_partitioning.swap(partitioning);
-             data_out.add_data_vector (static_partitioning, "partitioning");
-         }
-      }
-  }
-}
-
-
-template <int dim, typename VECTOR>
-void OutputProcessor<dim,VECTOR>::add_data_vector(const VECTOR &data_vector, 
-                                               const std::string &desc)
-{
-  AssertThrow(initialized, ExcNotInitialized());
-  deallog.push("AddingData");
-  std::vector<std::string> dd = Utilities::split_string_list(desc);
-  if(data_out.default_suffix() != "") {
-    if (dd.size() ==1 ) 
-      data_out.add_data_vector (data_vector, desc);
-    else 
-      data_out.add_data_vector (data_vector, dd);
-    deallog << "Added data: " << desc << std::endl;
-  }
-  deallog.pop();
-}
-
-
-template <int dim, typename VECTOR>
-void OutputProcessor<dim,VECTOR>::write_data_and_clear() {
-  AssertThrow(initialized, ExcNotInitialized());
-  AssertThrow(output_file, ExcIO());
-  deallog.push("WritingData");
-  if(data_out.default_suffix() != "") {
-      data_out.build_patches();
-    
-      if(n_mpi_processes > 1) {
-          data_out.write_deal_II_intermediate(output_file);
-      } else {
-          data_out.write(output_file);
-      }
-      deallog << "Wrote output file." << std::endl;
-      data_out.clear();
-      output_file.close();
-      deallog << "Reset output." << std::endl;
-  }
-  deallog.pop();
-}
-
-template <int dim, typename VECTOR>
-void OutputProcessor<dim,VECTOR>::dump_vector (const VECTOR &, 
-                                            const std::string & )
-{
-  Assert(false, ExcNotImplemented());
-  // Only specializations exist.
-}
diff --git a/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h b/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h
deleted file mode 100644 (file)
index 6a93c0b..0000000
+++ /dev/null
@@ -1,193 +0,0 @@
-/*
- * Immersed Boundary Problem:
- * 
- * Header Files
- *
- * Author: 
- * Luca Heltai <heltai@dimat.unipv.it>
- * =============
- * License: GPL.
- * =============
- * $Id: ibm_rhs.h,v 1.15 2005/04/05 14:46:49 luca Exp $
- *
- */
-#ifndef PARSED_SYMMETRIC_TENSOR_FUNCTION_H
-#define PARSED_SYMMETRIC_TENSOR_FUNCTION_H
-
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/function_parser.h>
-#include <deal.II/base/point.h>
-#include <deal.II/base/tensor.h>
-#include <deal.II/base/symmetric_tensor.h>
-
-#include <deal.II/lac/vector.h>
-
-using namespace dealii;
-
-template <int rank, int dim>
-class ParsedSymmetricTensorFunction
-{
- public:
-  /** Constructor. It initializes the function parser.*/
-  ParsedSymmetricTensorFunction();
-  
-  /** Declare parameters needed by this class. */
-  static void declare_parameters(ParameterHandler &prm);
-
-  /** Parse parameters needed by this class. */
-  void parse_parameters(ParameterHandler &prm);
-  
-  /** Get ONE value at the given point. Evaluate the parsed function
-      at the given point and return the symmetric tensor. */
-  const SymmetricTensor<rank, dim> & operator() (const Point<dim> &p) const;
-  
-  /** Get a vector of the values at a number of specified points */
-  void value_list(const std::vector< Point<dim> > &points,
-                 std::vector< SymmetricTensor<rank, dim> > &values);
-    
-  /** Get the tensor currently stored. */
-  inline const SymmetricTensor<rank, dim> & operator() () const {
-    return t;
-  };
-  
-  /** Get the size of the equivalent triangular matrix. If rank = 2,
-      this number is equal to dim, otherwise it is equal to dim*(dim+1)/2. */
-  static unsigned int get_dim_triangular();
-  
-  /** Set internal time of the function. */
-  void set_time(double);
- private:
-  FunctionParser<dim> f;
-  mutable SymmetricTensor<rank, dim> t;
-};
-
-
-template<int dim>
-inline 
-Tensor<2,dim> operator*(const SymmetricTensor<4,dim> & C, 
-                       const Tensor<2,dim> &gradv) {
-  Tensor<2,dim> T;
-  for(unsigned int i=0; i<dim; ++i)
-    for(unsigned int j=0; j<dim; ++j)
-      for(unsigned int k=0; k<dim; ++k)
-       for(unsigned int l=0; l<dim; ++l)
-         T[i][j] += C[i][j][k][l] * gradv[k][l];
-  return T;
-}
-
-
-template<int dim>
-inline 
-double operator*(const Tensor<2,dim> &gradv,
-                const SymmetricTensor<2,dim> & C) {
-  double T = 0;
-  for(unsigned int i=0; i<dim; ++i)
-    for(unsigned int j=0; j<dim; ++j)
-      T += C[i][j] * gradv[i][j];
-  return T;
-}
-
-
-template<int dim>
-inline 
-double double_contract(const SymmetricTensor<2,dim> & C,
-                      const Tensor<2,dim> &gradv) {
-  double T = 0;
-  for(unsigned int i=0; i<dim; ++i)
-    for(unsigned int j=0; j<dim; ++j)
-      T += C[i][j] * gradv[i][j];
-  return T;
-}
-
-
-
-template<int dim>
-inline 
-Tensor<1,dim> operator*(const SymmetricTensor<2,dim> & C,
-                       const Tensor<1,dim> &gradt) {
-  Tensor<1,dim> T;
-  for(unsigned int j=0; j<dim; ++j)
-    for(unsigned int i=0; i<dim; ++i)
-      T[j] += C[j][i] * gradt[i];
-  return T;
-}
-
-
-/* template<int dim> */
-/* inline  */
-/* double double_contract(const Tensor<2,dim> &gradv, */
-/*                    const Tensor<2,dim> & C) { */
-/*   double T = 0; */
-/*   for(unsigned int i=0; i<dim; ++i) */
-/*     for(unsigned int j=0; j<dim; ++j) */
-/*       T += C[i][j] * gradv[i][j]; */
-/*   return T; */
-/* } */
-
-
-
-/* template<int dim> */
-/* inline  */
-/* Tensor<2,dim> operator*(const SymmetricTensor<4,dim> & C,  */
-/*                     const Tensor<2,dim> &gradv) { */
-/*   Tensor<2,dim> T; */
-/*   for(unsigned int i=0; i<dim; ++i) */
-/*     for(unsigned int j=0; j<dim; ++j) */
-/*       for(unsigned int k=0; k<dim; ++k) */
-/*     for(unsigned int l=0; l<dim; ++l) */
-/*       T[i][j] += C[i][j][k][l] * gradv[k][l]; */
-/*   return T; */
-/* } */
-
-
-/* template<int dim> */
-/* inline  */
-/* double operator*(const Tensor<2,dim> &gradv, */
-/*              const SymmetricTensor<2,dim> & C) { */
-/*   double T = 0; */
-/*   for(unsigned int i=0; i<dim; ++i) */
-/*     for(unsigned int j=0; j<dim; ++j) */
-/*       T += C[i][j] * gradv[i][j]; */
-/*   return T; */
-/* } */
-
-
-/* template<int dim> */
-/* inline  */
-/* double double_contract(const SymmetricTensor<2,dim> & C, */
-/*                    const Tensor<2,dim> &gradv) { */
-/*   double T = 0; */
-/*   for(unsigned int i=0; i<dim; ++i) */
-/*     for(unsigned int j=0; j<dim; ++j) */
-/*       T += C[i][j] * gradv[i][j]; */
-/*   return T; */
-/* } */
-
-
-
-/* template<int dim> */
-/* inline  */
-/* Tensor<1,dim> operator*(const SymmetricTensor<2,dim> & C, */
-/*                     const Tensor<1,dim> &gradt) { */
-/*   Tensor<1,dim> T; */
-/*   for(unsigned int j=0; j<dim; ++j) */
-/*     for(unsigned int i=0; i<dim; ++i) */
-/*       T[j] += C[j][i] * gradt[i]; */
-/*   return T; */
-/* } */
-
-
-template<int dim>
-inline 
-double double_contract(const Tensor<2,dim> &gradv,
-                      const Tensor<2,dim> & C) {
-  double T = 0;
-  for(unsigned int i=0; i<dim; ++i)
-    for(unsigned int j=0; j<dim; ++j)
-      T += C[i][j] * gradv[i][j];
-  return T;
-}
-
-
-#endif
diff --git a/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.templates.h b/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.templates.h
deleted file mode 100644 (file)
index 3303319..0000000
+++ /dev/null
@@ -1,135 +0,0 @@
-/*
- * Immersed Boundary Problem:
- * 
- * Header Files
- *
- * Author: 
- * Luca Heltai <heltai@dimat.unipv.it>
- * =============
- * License: GPL.
- * =============
- * $Id: ibm_rhs.cc,v 1.34 2005/04/05 14:46:49 luca Exp $
- *
- */
-#include "../include/parsed_symmetric_tensor_function.h"
-#include <deal.II/base/utilities.h>
-
-template <int rank, int dim>
-  void ParsedSymmetricTensorFunction<rank, dim>::value_list(const std::vector< Point<dim> > &points,
-                                                           std::vector< SymmetricTensor<rank, dim> > &values)
-{
-
-  for (unsigned int n=0; n<points.size(); ++n) 
-    values[n] = this->operator()(points[n]);
-      
-}
-
-template <int rank, int dim>
-  unsigned int ParsedSymmetricTensorFunction<rank, dim>::get_dim_triangular() 
-{
-  if (rank == 2) 
-    return dim;
-  else if (rank == 4)
-    return (dim*(dim+1)/2);
-  else {
-    AssertThrow(false, ExcInternalError());
-    return 0;
-  }
-}
-
-
-template <int rank, int dim>
-void ParsedSymmetricTensorFunction<rank, dim>::set_time (double t) {
-  f.set_time(t);
-}
-
-
-template <int rank, int dim>
-ParsedSymmetricTensorFunction<rank, dim>::ParsedSymmetricTensorFunction () :
-  // Number of vector functions.
-  f(get_dim_triangular()*(get_dim_triangular()+1)/2)
-{
-  Assert(rank == 2 || rank == 4,
-        ExcMessage("Rank has to be even, either 2 or 4."));
-}
-
-template <int rank, int dim>
-void ParsedSymmetricTensorFunction<rank, dim>::declare_parameters(ParameterHandler &prm) 
-{
-  std::string vnames;
-  switch (dim) {
-  case 1:
-    vnames = "x,t";
-    break;
-  case 2:
-    vnames = "x,y,t";
-    break;
-  case 3:
-    vnames = "x,y,z,t";
-    break;
-  default:
-    AssertThrow(false, ExcNotImplemented());
-    break;
-  }
-  prm.declare_entry("Variable names", vnames, Patterns::Anything(), 
-                   "The name of the variables as they will be used in the function, separated by ','.");
-  // The expressions of the function
-  std::vector<std::string> expr(get_dim_triangular());
-  for(unsigned int i=0; i<get_dim_triangular(); ++i) {
-    for(unsigned int j=i; j<get_dim_triangular(); ++j)
-      expr[i] += (i==j) ? "1" : "; 0";
-    char tmp[100];
-    sprintf(tmp, "Row %d", i+1);
-    // Now we have the function expressions
-    prm.declare_entry(tmp, expr[i], Patterns::Anything(), (i == 0)  ?
-                     "Separate different components expressions by ';' "
-                     "as ',' is used internally by the function parser." 
-                     : "");
-  }
-  prm.declare_entry("Function constants", "", Patterns::Anything(),
-                   "Any constant used inside the functions which is not a variable name.");
-}
-
-template <int rank, int dim>
-  void ParsedSymmetricTensorFunction<rank, dim>::parse_parameters(ParameterHandler &prm) 
-{
-  std::string vnames = prm.get("Variable names");
-  std::string expr;
-  for(unsigned int i=0; i < get_dim_triangular(); ++i) {
-    char tmp[100];
-    sprintf(tmp, "Row %d", i+1);
-    expr += prm.get(tmp);
-    if(i+1 < get_dim_triangular())
-      expr += "; "; 
-  }
-  std::string constants_list = prm.get("Function constants");
-  
-  std::vector<std::string> const_list = 
-    Utilities::split_string_list(constants_list, ',');
-  std::map<std::string, double> constants;
-  for(unsigned int i = 0; i < const_list.size(); ++i) {
-    std::vector<std::string> this_c = 
-      Utilities::split_string_list(const_list[i], '=');
-    AssertThrow(this_c.size() == 2, ExcMessage("Invalid format"));
-    double tmp;
-    AssertThrow( sscanf(this_c[1].c_str(), "%lf", &tmp), ExcMessage("Double number?"));
-    constants[this_c[0]] = tmp;
-  }
-  
-  constants["pi"] = M_PI;
-  constants["Pi"] = M_PI;
-  
-  unsigned int nn = (Utilities::split_string_list(vnames)).size();
-  switch (nn) {
-  case dim:
-    // Time independent function
-    f.initialize(vnames, expr, constants); 
-    break;
-  case dim+1:
-    // Time dependent function
-    f.initialize(vnames, expr, constants, true);
-    break;
-  default:
-    AssertThrow(false, ExcMessage("Not the correct size. Check your code."));
-  }
-}
diff --git a/deal.II/examples/step-41/include/utilities.h b/deal.II/examples/step-41/include/utilities.h
deleted file mode 100644 (file)
index e1cdf83..0000000
+++ /dev/null
@@ -1,62 +0,0 @@
-#ifndef UTILITIES_HLT
-#define UTILITIES_HLT
-
-#include <deal.II/base/utilities.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/lac/vector_memory.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-
-using namespace dealii;
-
-template <class Matrix>
-class InverseMatrix : public Subscriptor
-{
- public:
-  InverseMatrix (const Matrix &m);
-  
-  void vmult (Vector<double>       &dst,
-             const Vector<double> &src) const;
-  
-  void Tvmult (Vector<double>       &dst,
-              const Vector<double> &src) const;
- private:
-  const SmartPointer<const Matrix> matrix;
-  mutable GrowingVectorMemory<> vector_memory;    
-};
-  
-class BBt : public Subscriptor
-{
- public:
-  void reinit (const BlockSparseMatrix<double> &A);
-  
-  
-  void vmult (Vector<double>       &dst,
-              const Vector<double> &src) const;
-
-  void Tvmult (Vector<double>       &dst,
-              const Vector<double> &src) const;
-
- private:
-  SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-
-  mutable BlockVector<double> u, v;
-
-  unsigned int dim;
-};
-
-template <typename TYPE>
-void smart_delete (SmartPointer<TYPE> &sp) {
-  if(sp) {
-    TYPE * p = sp;
-    sp = 0;
-    delete p;
-  }
-}
-
-#endif
diff --git a/deal.II/examples/step-41/include/vector_space.h b/deal.II/examples/step-41/include/vector_space.h
deleted file mode 100644 (file)
index 1d200f7..0000000
+++ /dev/null
@@ -1,325 +0,0 @@
-#ifndef VECTOR_SPACE_H
-#define VECTOR_SPACE_H
-
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/grid/persistent_tria.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/fe/mapping.h>
-#include <deal.II/dofs/dof_constraints.h>
-#include <deal.II/numerics/solution_transfer.h>
-
-using namespace dealii;
-
-/** 
-  VectorSpace object. Anything related to the definition of the
-  Hilbert space that makes up the pde problem lives in this class.
-*/
-template <int dim>
-class VectorSpace : public Subscriptor
-{
-  public:
-  
-  /** Empty constructor. */
-  VectorSpace ();
-  /** Full constructor. It initializes this object by reading a
-      parameter file and setting pointers to the given
-      triangulation. */
-  VectorSpace (ParameterHandler &prm, 
-              Triangulation<dim> &dd,
-              const unsigned int n_mpi_processes=1,
-              const unsigned int this_mpi_process=0);
-
-  /** Destructor. Clears the created pointers. */
-  ~VectorSpace ();
-  
-  /** Reinit. The state of this object after calling reinit, is
-      the same as after the full constructor is used.*/
-  void reinit(ParameterHandler &prm, 
-             Triangulation<dim> &dd,
-             const unsigned int n_mpi_processes=1,
-             const unsigned int this_mpi_process=0,
-             const std::string space_name="Vector Space Parameters");
-
-  /** Reset the internal status of the current triangulation, leaving
-      untouched all the rest. The status of this object is as if the
-      current triangulation and dof handler had just been
-      constructed. Nothing is done to the saved triangulation, nor to
-      the saved dofhandler. */
-  void reinit();
-  
-  /** Parse the parameter file. */
-  void parse_parameters(ParameterHandler &prm, const std::string space_name);
-
-  /** Generate entries in the given parameter file. */
-  static void declare_parameters(ParameterHandler &prm, const std::string space_name="Vector Space Parameters");
-
-  /** Initialize the mesh. */
-  void initialize_mesh();
-
-  /** Refine the mesh globally. */
-  void refine_globally();
-  
-  /** Flag the mesh according to the given error estimator and to the
-      refinement strategy defined in the paratmeters. */
-  void flag_for_refinement(Vector<float> &error);
-  
-  /** Save current vector space. */
-  void save_step();
-  
-  /** Save current vector space, and reset triangulation to be pristine. */
-  void save_step_and_reinit();
-
-  /** Transfer solution from old to new grid. */
-  template <typename VECTOR>
-  void transfer(VECTOR &dst, const VECTOR &src);
-
-  /** Reorder the degrees of freedom according to the
-      parameters. This is done for the level specified on the
-      argument. If level is -1, then the active dofs are reordered. */
-  void reorder(int level = -1,
-              std::vector<unsigned int> target_comps=std::vector<unsigned int>());
-
-  /** Redistribute degrees of freedom. The optional arguments groups
-      together the components. It is useful for vector valued finite elements, when one wants to sort together*/
-  void redistribute_dofs(std::vector<unsigned int> 
-                        target_components=std::vector<unsigned int>());
-      
-  /** Compute maximum and minimum diameter of the mesh cells. */
-  void measure_mesh(); 
-      
-  /** Interpolate Boundary Conditions. Generates the boundary
-      conditions for the problem. This will be used after the assembly
-      procedures. */
-  void interpolate_dirichlet_bc(const Function<dim> & f,  
-                               std::map<unsigned int, double> & bvalues);
-   
-  /** Return reference to finite element.*/
-  inline FiniteElement<dim,dim> & get_fe() {
-    return *fe;
-  }
-
-  /** Return reference to current dof handler.*/
-  inline MGDoFHandler<dim,dim> & get_dh() {
-    return *dh;
-  }
-
-  /** Return reference to previous dof handler.*/
-  inline MGDoFHandler<dim,dim> & get_last_dh() {
-      return *last_dh;
-  }
-
-  /** Return reference to current triangulation..*/
-  inline Triangulation<dim,dim> & get_tria() {
-    return *tria;
-  }
-
-  /** Return reference to current triangulation..*/
-  inline Triangulation<dim> & get_last_tria() {
-      return *last_tria;
-  }
-
-  /** Return reference to coarse triangulation..*/
-  inline Triangulation<dim> & get_coarse_tria() {
-    return *coarse;
-  }
-
-  /** Return reference to hanging node constraints..*/
-  inline ConstraintMatrix & get_hang() {
-    return hang;
-  }
-      
-  /** Return reference to mapping.*/
-  inline Mapping<dim> & get_mapping() {
-    return *mapping;
-  }
-
-      
-  /** Return constant reference to finite element.*/
-  inline const FiniteElement<dim,dim> & get_fe() const {
-    return *fe;
-  }
-
-  /** Return constant reference to current dof handler.*/
-  inline const MGDoFHandler<dim,dim> & get_dh() const {
-    return *dh;
-  }
-
-  /** Return constant reference to previous dof handler.*/
-  inline const MGDoFHandler<dim,dim> & get_last_dh() const {
-      return *last_dh;
-  }
-
-  /** Return reference to current triangulation..*/
-  inline const Triangulation<dim,dim> & get_tria() const {
-    return *tria;
-  }
-
-  /** Return reference to coarse triangulation..*/
-  inline const Triangulation<dim,dim> & get_coarse_tria() const {
-    return *coarse;
-  }
-
-  /** Return reference to current triangulation..*/
-  inline const Triangulation<dim,dim> & get_last_tria() const {
-      return *last_tria;
-  }
-
-  /** Return reference to hanging node constraints..*/
-  inline const ConstraintMatrix & get_hang() const {
-    return hang;
-  }
-      
-  /** Return reference to mapping.*/
-  inline const Mapping<dim,dim> & get_mapping() const {
-    return *mapping;
-  }
-  
-  /** Number of dofs. */
-  inline unsigned int n_dofs() {
-    return dh->n_dofs();
-  }
-  
-  /** Number of dofs per process. */
-  inline const std::vector<unsigned int> & n_dofs_pp() const {
-    return local_dofs_per_process;
-  }
-  
-  /** Number of dofs per block. */
-  inline const std::vector<unsigned int> & n_dofs_pb() const {
-    return dofs_per_block;
-  }
-
-  /** Number of dofs in this process. */
-  inline  unsigned int n_ldofs() const {
-    return n_local_dofs;
-  }
-  
-  /** Number of dofs per block per process.  */
-  inline const std::vector<std::vector<unsigned int> > & n_dofs_pb_pp() const {
-    return dofs_per_block_per_process;
-  }
-
-  /** Number of dofs per process per block.  */
-  inline const  std::vector<std::vector<unsigned int> > & n_dofs_pp_pb() const {
-    return dofs_per_process_per_block;
-  }
-
-  /** Number of blocks. */
-  inline  unsigned int n_blocks() const {
-    return number_of_blocks;
-  }
-  /** Bool that checks if local refiment is enabled.*/
-  bool enable_local_refinement;
-  
-  /** Initial global refinement. */
-  unsigned int global_refinement;
-
-  /** The size of this mesh. */
-  double h_max;
-      
-  /** The size of this mesh. */
-  double h_min;
-      
-  /** Dirichlet Boundary Indicators. */
-  std::map<char, std::vector<bool> > dirichlet_bc;
-      
-  /** Neumann Boundary Indicators. */
-  std::map<char, std::vector<bool> > neumann_bc;
-      
-  /** Other Boundary Indicators. */
-  std::map<char, std::vector<bool> > other_bc;
-private:
-  
-  /** For internal use. Counts the dofs per block, block and what not. */
-  void count_dofs(std::vector<unsigned int> target_components);
-  
-  std::string fe_name;
-  std::vector<std::string> ordering;
-  unsigned int map_type;
-
-  std::vector<std::string> d_bc;
-  std::vector<std::string> n_bc;
-  std::vector<std::string> o_bc;
-      
-  /** Pointer to the finite element used. */
-  SmartPointer<FiniteElement<dim,dim> > fe;
-
-  /** Pointer to the dofhandler used */
-  SmartPointer<MGDoFHandler<dim,dim> > dh;
-
-  /** Pointer to the last dofhandler used */
-  SmartPointer<MGDoFHandler<dim,dim> > last_dh;
-
-  /** Pointer to a pristine coarse triangulation. */
-  SmartPointer<Triangulation<dim,dim> > coarse;
-
-  /** Pointer to the current triangulation used */
-  SmartPointer<Triangulation<dim,dim> > tria;
-
-  /** Pointer to the previous triangulation used. */
-  SmartPointer<Triangulation<dim,dim> > last_tria;
-
-  /** Finite Element Mapping. Various mappings are supported. If the 
-      mapping parameter is 0, then cartesian mapping is used. Else Qn, with 
-      n the mapping parameter. */
-  SmartPointer<Mapping<dim,dim> > mapping;
-
-  /** Constraint Matrix. */
-  ConstraintMatrix hang;
-
-  public:
-      
-  /** Number of processes. */
-  unsigned int n_mpi_processes;
-
-  /** The id of this process. */
-  unsigned int this_mpi_process;
-
-  private:
-  
-  /** Local size of vectors. */
-  std::vector<unsigned int> local_dofs_per_process;
-
-  /** Local size of vectors. */
-  std::vector<unsigned int> dofs_per_block;
-
-  /** Local sizes of vectors divided by blocks. */
-  std::vector<std::vector<unsigned int> > dofs_per_block_per_process;
-
-  /** Local sizes of vectors divided by blocks. */
-  std::vector<std::vector<unsigned int> > dofs_per_process_per_block;
-
-  /** Number of local dofs. */
-  unsigned int n_local_dofs;
-  
-  /** Number of Blocks. */
-  unsigned int number_of_blocks;
-
-  /** Number of cells. */
-  unsigned int number_of_cells;
-  
-  /** Distortion coefficient. */
-  double distortion;
-  
-  /** Refinement strategy. */
-  std::string refinement_strategy;
-
- public:
-
-  /** Bottom fraction of refinement. */
-  double bottom_fraction;
-
-  /** Top fraction of refinement. */
-  double top_fraction;
-
- private:
-
-  /** Maximum number of allowed cells. */
-  unsigned int max_cells;
-    
-  /** The wind direction, in case we order the mesh upwind. */
-  Point<dim> wind;
-};
-#endif
diff --git a/deal.II/examples/step-41/include/vector_space.templates.h b/deal.II/examples/step-41/include/vector_space.templates.h
deleted file mode 100644 (file)
index 1a10fb3..0000000
+++ /dev/null
@@ -1,583 +0,0 @@
-#include "../include/vector_space.h"
-#include "../include/utilities.h"
-#include <deal.II/base/logstream.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_in.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/fe/mapping_cartesian.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/dofs/dof_renumbering.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/numerics/solution_transfer.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/grid_refinement.h>
-
-using namespace std;
-
-template <int dim> 
-VectorSpace<dim>::VectorSpace() :
-    fe(0, "Vector Space FE"),
-    dh(0, "Vector Space This DH"),
-    last_dh(0, "Vector Space Last DH"),
-    coarse(0, "Vector Space Coarse Grid"),
-    tria(0, "Vector Space This Grid"),
-    last_tria(0, "Vector Space Last Grid"),
-    mapping(0, "Vetor Space Mapping")
-{}
-
-template <int dim>
-VectorSpace<dim>::VectorSpace(ParameterHandler &prm, 
-                             Triangulation<dim> &dd,
-                             const unsigned int n_mpi,
-                             const unsigned int this_mpi) :
-    fe(0, "Vector Space FE"),
-    dh(0, "Vector Space This DH"),
-    last_dh(0, "Vector Space Last DH"),
-    coarse(0, "Vector Space Coarse Grid"),
-    tria(0, "Vector Space This Grid"),
-    last_tria(0, "Vector Space Last Grid"),
-    mapping(0, "Vetor Space Mapping"),
-    n_mpi_processes(n_mpi),
-    this_mpi_process(this_mpi)
-{
-  reinit(prm, dd, n_mpi, this_mpi);
-}
-
-template <int dim> 
-VectorSpace<dim>::~VectorSpace()
-{
-  // Get rid of the dof handler object.
-  smart_delete(dh);
-  smart_delete(last_dh);
-  smart_delete(fe);
-  coarse = 0;
-  smart_delete(tria);
-  smart_delete(last_tria);
-  smart_delete(mapping);
-}
-
-template <int dim>
-void VectorSpace<dim>::save_step() {
-  deallog.push("Save");
-  // Get rid of old dh and triangulation
-  last_dh->clear();
-  last_tria->clear();
-  
-  deallog << "Saving triangulation and dofs..." << std::endl;
-  // Copy old triangulation (which is the one just used)
-  last_tria->copy_triangulation(*tria);
-  
-  // swap back and forth, so that we have the same dofs...
-  dh.swap(last_dh);
-  redistribute_dofs();
-  dh.swap(last_dh);
-  deallog.pop();
-}
-
-
-template <int dim>
-void VectorSpace<dim>::save_step_and_reinit() {
-  deallog.push("Save");
-  deallog << "Swapping tria and dh with old_tria and old_dh." << std::endl;
-  // Swap dofhandler.
-  last_dh.swap(dh);
-  // Swap triangulation.
-  last_tria.swap(tria);
-  // Now call reinit
-  reinit();
-  deallog.pop();
-}
-
-
-template <int dim>
-void VectorSpace<dim>::refine_globally() 
-{
-    deallog.push("GlobalRefinement");
-    tria->refine_global (1);
-    redistribute_dofs();
-    deallog.pop();
-}
-
-template <int dim>
-void VectorSpace<dim>::flag_for_refinement(Vector<float> &error) 
-{
-    deallog.push("Refining");
-    deallog << "Strategy: " << refinement_strategy << std::endl;
-    if(refinement_strategy == "global") {
-      tria->set_all_refine_flags();
-    } else if(refinement_strategy == "fixed_number") {
-      if(max_cells > 0)
-       GridRefinement::refine_and_coarsen_fixed_number(*tria, error, 
-                                                       top_fraction, bottom_fraction, max_cells);
-      else
-       GridRefinement::refine_and_coarsen_fixed_number(*tria, error, 
-                                                       top_fraction, bottom_fraction);
-    } else if(refinement_strategy == "fixed_fraction") {
-      if(max_cells > 0)
-       GridRefinement::refine_and_coarsen_fixed_fraction(*tria, error, 
-                                                         top_fraction, bottom_fraction, max_cells);
-      else
-       GridRefinement::refine_and_coarsen_fixed_fraction(*tria, error, 
-                                                         top_fraction, bottom_fraction);
-    } else if(refinement_strategy == "optimize") {
-      GridRefinement::refine_and_coarsen_optimize(*tria, error);
-    } else {
-      Assert(false, ExcInternalError());
-    }
-    deallog.pop();
-}
-
-
-
-template <int dim>
-void VectorSpace<dim>::reorder(int level,
-                              std::vector<unsigned int> target_comps) {
-
-  if(level<0) level = dh->get_tria().n_levels()-1;
-
-  // Renumber by subdomain association only if needed
-  if(n_mpi_processes > 1) {
-    deallog << "Renumbering subdomain wise." << std::endl;
-    DoFRenumbering::subdomain_wise(static_cast<DoFHandler<dim>&>(*dh));
-  }
-  
-  for(unsigned int i=0; i<ordering.size(); ++i) {
-    if(ordering[i] == "cuth") {
-      deallog << "Renumbering with Cuthill McKee algorithm." << std::endl;
-      DoFRenumbering::Cuthill_McKee(static_cast<DoFHandler<dim>&>(*dh));
-    } else if(ordering[i] == "comp") {
-      deallog << "Renumbering component wise." << std::endl;
-      DoFRenumbering::component_wise(static_cast<DoFHandler<dim>&>(*dh), target_comps);
-    } else if(ordering[i] == "upwind") {
-      DoFRenumbering::downstream(static_cast<DoFHandler<dim>&>(*dh), wind);
-    } else if(ordering[i] == "none") {}
-    else {
-      AssertThrow(false, ExcMessage("This Reordering not implemented"));
-    }
-  }
-    
-}
-
-template <int dim>
-void VectorSpace<dim>::measure_mesh() {
-  typename Triangulation<dim>::active_cell_iterator cell, endc;
-  endc = tria->end();
-  h_max = 0;
-  h_min = 1000;
-  for(cell = tria->begin_active(); cell != endc; ++cell) {
-    h_max = std::max(h_max, cell->diameter());
-    h_min = std::min(h_min, cell->diameter());
-  }
-  deallog << "Max diameter of a cell: " << h_max << std::endl
-         << "Min diameter of a cell: " << h_min << std::endl;
-}
-
-template <int dim>
-void VectorSpace<dim>::reinit(ParameterHandler &prm, 
-                             Triangulation<dim> &dd,
-                             const unsigned int n_mpi,
-                             const unsigned int this_mpi,
-                             const std::string space_name) 
-{
-  deallog.push("VECTORSPACE");
-  n_mpi_processes = n_mpi;
-  this_mpi_process = this_mpi;
-  
-  smart_delete(dh);
-  smart_delete(last_dh);
-  smart_delete(fe);
-  smart_delete(mapping);
-  smart_delete(coarse);
-  smart_delete(last_tria);
-
-  // Parse Parameters
-  parse_parameters(prm, space_name);
-
-  // generate Finite element
-  FiniteElement<dim> * fe_p=0;
-    
-  try {
-      fe_p = FETools::get_fe_from_name<dim>(fe_name);
-  } catch (...) {
-
-    // Did not recognize the finite element to use. Try your own
-    // one...
-    // fe_p = new FESystem<dim>(FE_BGH<dim>(2),1, FE_DGP<dim>(1), 1);
-    throw;
-  }
-  fe = fe_p;
-  deallog << "Finite Element Space: " << fe->get_name() << endl;
-
-  if(map_type == 0) {
-    mapping= new MappingCartesian<dim>();
-  } else {
-    mapping = new MappingQ<dim>(map_type);
-  }
-    
-  // Store a pointer to the triangulation
-  coarse = &dd;
-
-  // Now generate an empty triangulation.
-  last_tria = new Triangulation<dim>();
-  
-  // And a copy of the pristine version of the triangulation.
-  tria = new Triangulation<dim>();
-  tria->copy_triangulation(*coarse);
-
-  // And we initialize the current triangulation.
-  initialize_mesh();
-
-  // Generate new DoFHandlers
-  dh = new MGDoFHandler<dim>(*tria);
-  last_dh = new MGDoFHandler<dim>(*last_tria);
-  
-  // Now generates the boundary maps
-  std::vector<std::vector<std::string> > bcs;
-  bcs.push_back(d_bc);
-  bcs.push_back(n_bc);
-  bcs.push_back(o_bc);
-
-  std::vector<std::map<char, std::vector<bool> > > bc_maps(3);
-    
-  for(unsigned int bcnumber=0; bcnumber<3; ++bcnumber) {
-    // Get some aliases for the bcs
-    std::vector<std::string> & bc = bcs[bcnumber];
-    std::map<char, std::vector<bool> > & bc_map = bc_maps[bcnumber];
-      
-    if(bc.size()) 
-      for(unsigned int i=0; i < bc.size(); ++i) {
-       std::vector<std::string> id_and_comps = 
-         Utilities::split_string_list(bc[i], ':');
-       AssertThrow(id_and_comps.size() == 2,
-                   ExcMessage("Wrong Format for boundary indicator map."));
-       std::vector<int> ids = 
-         Utilities::string_to_int(Utilities::split_string_list(id_and_comps[0]));
-       std::vector<int> comps = 
-         Utilities::string_to_int(Utilities::split_string_list(id_and_comps[1]));
-         
-       unsigned int n_c = fe->n_components();
-       std::vector<bool> filter(n_c, false);
-       // Now check that the components make sense
-       for(unsigned int i=0; i<comps.size(); ++i) {
-         AssertThrow((unsigned int) comps[i] < n_c,
-                     ExcIndexRange(comps[i], 0, n_c));
-         filter[ comps[i] ] = true;
-       }
-       // And now save these components. Merge them if the map is
-       // already stored.
-       for(unsigned int i = 0; i<ids.size(); ++i) {
-         // For each of the id, save the map just generated
-         if(bc_map.find(ids[i]) != bc_map.end()) {
-           // For each of the id, add the components that appear in
-           // the map just generated.
-           Assert(bc_map[ ids[i] ].size() == (unsigned int) n_c,
-                  ExcDimensionMismatch(bc_map[ ids[i] ].size(), n_c));
-           for(unsigned int j=0; j<n_c; ++j) {
-             bc_map[ ids[i] ][j] = filter[j];
-           }
-         } else {
-           bc_map[ids[i]] = filter;
-         }
-       }
-      }
-  }
-  dirichlet_bc = bc_maps[0];
-  neumann_bc = bc_maps[1];
-  other_bc = bc_maps[2];
-
-  deallog.pop();
-}
-
-
-
-template <int dim>
-void VectorSpace<dim>::reinit() {
-  // In this case, we only get rid of the current dofhandler
-  // and triangulation, assuming they have already been saved.
-  dh->clear();
-  tria->clear();
-  // And we restore the initial copy of the triangulation.
-  tria->copy_triangulation(*coarse);
-}
-
-template <int dim>
-void VectorSpace<dim>::initialize_mesh() 
-{
-  deallog.push("InitializeMesh");
-  tria->refine_global (global_refinement);
-  deallog << "Active Cells: "
-         << tria->n_active_cells()
-         << endl;
-  deallog.pop();
-}
-
-template <int dim>
-void VectorSpace<dim>::redistribute_dofs(std::vector<unsigned int> target_components ) {
-    // Repartition the triangulation
-    if(n_mpi_processes > 1) 
-      GridTools::partition_triangulation (n_mpi_processes, *tria);
-    
-    // Distort the mesh if necessary
-    if(distortion > 0.) tria->distort_random(distortion);
-    
-    // Measure it
-    measure_mesh();
-    
-    // And Distribute degrees of freedom
-    dh->distribute_dofs(*fe);
-    
-    // Now count what you just distributed
-    count_dofs(target_components);
-    
-    reorder(-1, target_components);
-    hang.clear();
-    DoFTools::make_hanging_node_constraints (*dh, hang);
-    hang.close();
-    deallog << "Number of constrained DOFS: " << hang.n_constraints() << std::endl;
-}
-
-template <int dim>
-void VectorSpace<dim>::declare_parameters(ParameterHandler &prm, const std::string space_name) 
-{
-
-   prm.enter_subsection(space_name);
-
-    prm.declare_entry ("Finite element space", "FE_Q(1)", 
-                      Patterns::Anything(),
-                      "The finite element space to use. For vector "
-                      "finite elements use the notation "
-                      "FESystem[FE_Q(2)^2-FE_DGP(1)] (e.g. Navier-Stokes). ");
-
-    prm.declare_entry ("Mapping degree", "1", Patterns::Integer(),
-                      "Degree of the mapping. If 0 is used, then a Cartesian mapping is assumed.");
-    prm.declare_entry ("Dof ordering", "cuth, comp", Patterns::Anything(),
-                      "Ordering of the degrees of freedom: none, comp, cuth, upwind.");
-    prm.declare_entry ("Wind direction", ".01, .01, 1", Patterns::Anything(),
-                      "Direction of the wind for upwind ordering of the mesh. ");
-
-    prm.declare_entry ("Dirichlet boundary map", "1:0", Patterns::Anything(),
-                      "Boundary indicator, followed by semicolomn and a list"
-                      " of components to which this boundary conditions apply. "
-                      "More boundary indicators can be separated by semicolumn. "
-                      "1:0,1,4 ; 2,4:0,2");
-    prm.declare_entry ("Neumann boundary map", "2:0", Patterns::Anything(),
-                      "Boundary indicators, followed by semicolomn and a list of "
-                      "components to which this boundary conditions apply. More "
-                      "boundary indicators can be separated by semicolumn. "
-                      "1:0,1,4 ; 2,4:0,2");
-    prm.declare_entry ("Other boundary map", "3:0", Patterns::Anything(),
-                      "Boundary indicator, followed by semicolomn and a list of "
-                      "components to which this boundary conditions apply. More "
-                      "boundary indicators can be separated by semicolumn. "
-                      "1:0,1,4 ; 2,4:0,2");
-    
-    prm.enter_subsection("Grid Parameters");
-    
-    prm.declare_entry ("Global refinement", "4", Patterns::Integer());
-    prm.declare_entry ("Distortion coefficient", "0", Patterns::Double(),
-                      "If this number is greater than zero, the mesh is distorted"
-                      " upon refinement in order to disrupt its structureness.");
-    
-    prm.declare_entry("Refinement strategy", 
-                     "fixed_number", Patterns::Selection("fixed_number|fixed_fraction|optimize|global"),
-                     "fixed_number: the Top/Bottom threshold fraction of cells are flagged for "
-                     "refinement/coarsening. "
-                     "fixed_fraction: the cells whose error is Top/Bottom fraction of the total "
-                     "are refined/coarsened. optmized: try to reach optimal error distribution, "
-                     "assuming error is divided by 4 upon refining. global: refine all cells.");
-    prm.declare_entry("Bottom fraction", ".3", Patterns::Double());
-    prm.declare_entry("Top fraction", ".3", Patterns::Double());
-    prm.declare_entry("Max number of cells", "0", Patterns::Integer(),
-                     "A number of zero means no limit. ");
-    prm.leave_subsection();
-  
-    
-    prm.leave_subsection();
-}
-
-    template <int dim>
-void VectorSpace<dim>::parse_parameters(ParameterHandler &prm, const std::string space_name) 
-{
-    prm.enter_subsection(space_name);
-
-    fe_name = prm.get ("Finite element space"); 
-    map_type = prm.get_integer("Mapping degree");
-    std::string all_ordering = prm.get ("Dof ordering");
-
-    d_bc =  Utilities::split_string_list(prm.get ("Dirichlet boundary map"), ';');
-    n_bc =  Utilities::split_string_list(prm.get ("Neumann boundary map"), ';');
-    o_bc =  Utilities::split_string_list(prm.get ("Other boundary map"), ';');
-    
-    std::vector<std::string> wind_str = 
-       Utilities::split_string_list(prm.get ("Wind direction") );
-    for(unsigned int i=0; (i<wind_str.size()) && (i<dim); ++i) 
-           sscanf(wind_str[i].c_str(), "%lf", &wind[i]);
-                                    
-    prm.enter_subsection("Grid Parameters");
-    
-    global_refinement          = prm.get_integer ("Global refinement");
-    distortion                 = prm.get_double("Distortion coefficient");
-    
-    refinement_strategy                = prm.get("Refinement strategy");
-    enable_local_refinement    = !(refinement_strategy == "global");
-    
-    bottom_fraction            = prm.get_double("Bottom fraction");
-    top_fraction               = prm.get_double("Top fraction");
-    max_cells                  = prm.get_integer("Max number of cells");
-    prm.leave_subsection();
-
-    prm.leave_subsection();
-
-    ordering = Utilities::split_string_list(all_ordering);
-}
-
-
-template <int dim>
-void VectorSpace<dim>::interpolate_dirichlet_bc(const Function<dim> & f,  
-                                               std::map<unsigned int, double> & bvalues) 
-{
-  deallog.push("DBC");
-  std::map<char, std::vector<bool> >::iterator 
-    dmap = dirichlet_bc.begin(),
-    dmapend = dirichlet_bc.end();
-  
-  unsigned int last_counted = 0;
-  for(; dmap!=dmapend; ++dmap) {
-    char id = dmap->first;
-    std::vector<bool> &filter = dmap->second;
-    
-        deallog << (int) id << " :";
-    const unsigned int n_components = fe->n_components();
-    for(unsigned int i=0; i < n_components; ++i) {
-      if(filter[i]) ;//deallog << i << ", ";
-    }
-    VectorTools::interpolate_boundary_values(*mapping, (const DoFHandler<dim>&) *dh, 
-                                            // Dirichlet boundary only...
-                                            id, f, bvalues, filter);
-    deallog << " #: " << bvalues.size() - last_counted << std::endl;
-    last_counted = bvalues.size();
-  }
-  deallog << "Total DBC dofs: " << bvalues.size() << std::endl;
-  deallog.pop();
-}
-
-
-
-template <int dim>
-void VectorSpace<dim>::count_dofs(std::vector<unsigned int> target_components) {
-  // Count dofs per processor
-  local_dofs_per_process.resize(n_mpi_processes);
-  
-  // Output dofs per processor
-  deallog << "Number of DOFS: " << dh->n_dofs() 
-         << ", (by partition:";
-  for (unsigned int p=0; p<n_mpi_processes; ++p) {
-    unsigned int nl = DoFTools::
-      count_dofs_with_subdomain_association (*dh, p);
-    local_dofs_per_process[p] = nl;
-    deallog << (p==0 ? ' ' : '+') << nl;
-  }    
-  deallog << ")" << std::endl;
-  
-  // Count number of blocks
-  number_of_blocks = ( target_components.size() ? 
-                      target_components[target_components.size()-1]+1 : 1 );
-  
-  deallog << "Number of blocks: " << number_of_blocks << std::endl;
-  
-  if( number_of_blocks > 1 ) {
-    // Count dofs per blocks.
-    dofs_per_block.resize( number_of_blocks );
-
-    DoFTools::count_dofs_per_component (*dh, dofs_per_block, false, target_components);
-    
-    while(dofs_per_block.size() > number_of_blocks) {
-      deallog << "Old Deal.II!!!!"<< std::endl;
-      dofs_per_block.pop_back();
-    }
-      
-    // Output the dofs per block
-    deallog << "Dofs by Block :";
-    for (unsigned int p=0; p<number_of_blocks; ++p) {
-      deallog << (p==0 ? ' ' : '+') << dofs_per_block[p];
-    }    
-    deallog << std::endl;
-    
-    // Now count dofs per block per processor
-    dofs_per_block_per_process.resize(number_of_blocks,
-                                     std::vector<unsigned int>(n_mpi_processes, 0));
-    dofs_per_process_per_block.resize(n_mpi_processes,
-                                     std::vector<unsigned int>(number_of_blocks, 0));
-
-    // We now have a vector that contains all subdomain ids for each dof
-    std::vector<unsigned int> subdomain_association(dh->n_dofs());
-    DoFTools::get_subdomain_association (*dh, subdomain_association);
-    
-    // Now cycle on each block to find out who belongs here
-    for(unsigned int block = 0; block < number_of_blocks; ++block) {
-      
-      // Select dofs belonging to this block
-      deallog << "Block " << block << ":";
-      std::vector<bool> mask(fe->n_components(), false);
-      for(unsigned int i=0; i<mask.size(); ++i) {
-       if(target_components[i] == block) mask[i] = true;
-       deallog << ( i==0 ? ' ' : ',' ) << mask[i];
-      }
-      deallog << ": ";
-      
-      std::vector<bool> this_block_dofs(dh->n_dofs());
-      DoFTools::extract_dofs(*dh, mask, this_block_dofs);
-      
-      // Reset counter
-      for(unsigned int j=0; j<n_mpi_processes; ++j) {
-       dofs_per_block_per_process[block][j] = 0;
-       dofs_per_process_per_block[j][block] = 0;
-      }
-      
-      // Count dofs on this block
-      for(unsigned int i=0; i< dh->n_dofs(); ++i)
-       if(this_block_dofs[i]) {
-         ++dofs_per_block_per_process[block][subdomain_association[i]];
-         ++dofs_per_process_per_block[subdomain_association[i]][block];
-       }
-      unsigned int check_pb=0;
-      unsigned int check_pp=0;
-      
-      // Output the dofs per block per processor
-      for(unsigned int j=0; j<n_mpi_processes; ++j) {
-       deallog << (j==0 ? ' ' : '+') << dofs_per_block_per_process[block][j];
-       check_pb += dofs_per_block_per_process[block][j];
-       check_pp += dofs_per_process_per_block[j][block];
-      }    
-      deallog << std::endl;
-      Assert(check_pb == dofs_per_block[block],
-            ExcDimensionMismatch(check_pb, dofs_per_block[block]));
-      Assert(check_pp == dofs_per_block[block],
-            ExcDimensionMismatch(check_pp, dofs_per_block[block]));
-    }
-  } else {
-    //
-    dofs_per_block.resize(1);
-    dofs_per_block[0] = dh->n_dofs();
-    // 
-    dofs_per_block_per_process.resize(1);
-    dofs_per_block_per_process[0] = local_dofs_per_process;
-    dofs_per_process_per_block.resize(n_mpi_processes);
-    for(unsigned int j=0; j<n_mpi_processes; ++j)
-      dofs_per_process_per_block[j] = 
-       std::vector<unsigned int>(1, local_dofs_per_process[j]);
-      
-  }
-  
-  // Save the number of cells.
-  number_of_cells = tria->n_active_cells();
-  deallog << "Number of CELLS: " << number_of_cells << endl;
-  
-  // Now save the number of local dofs...
-  n_local_dofs = local_dofs_per_process[this_mpi_process];
-}
diff --git a/deal.II/examples/step-41/mesh/square.inp b/deal.II/examples/step-41/mesh/square.inp
deleted file mode 100644 (file)
index 23f8127..0000000
+++ /dev/null
@@ -1,16 +0,0 @@
-# This file was generated by the deal.II library.
-# Date =  2006/1/6
-# Time =  18: 2:34
-#
-# For a description of the UCD format see the AVS Developer's guide.
-#
-4 5 0 0 0
-1  0 0 0
-2  1 0 0
-3  0 1 0
-4  1 1 0
-1 0 quad    1 2 4 3 
-2  1  line    1 2 
-3  4  line    1 3 
-4  2  line    2 4 
-5  3  line    3 4 
diff --git a/deal.II/examples/step-41/source/base.cc b/deal.II/examples/step-41/source/base.cc
deleted file mode 100644 (file)
index 9da975a..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-#include "../include/base.templates.h"
-
-template class Base<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/camclay.cc b/deal.II/examples/step-41/source/camclay.cc
deleted file mode 100644 (file)
index 506ae0c..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-#include "../include/camclay.templates.h"
-
-template class CamClay<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/domain.cc b/deal.II/examples/step-41/source/domain.cc
deleted file mode 100644 (file)
index fa3bc74..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-#include "../include/domain.templates.h"
-
-template class Domain<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/error_handler.cc b/deal.II/examples/step-41/source/error_handler.cc
deleted file mode 100644 (file)
index f4a7e5a..0000000
+++ /dev/null
@@ -1,11 +0,0 @@
-#include "../include/error_handler.templates.h"
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/petsc_vector.h>
-#include <deal.II/lac/vector.h>
-
-template class ErrorHandler<deal_II_dimension, Vector<double> >;
-template class ErrorHandler<deal_II_dimension, BlockVector<double> >;
-#ifdef CONTRIB_USE_PETSC
-template class ErrorHandler<deal_II_dimension, PETScWrappers::Vector >;
-#endif
-//template class ErrorHandler<deal_II_dimension, PETScWrappers::BlockVector >;
diff --git a/deal.II/examples/step-41/source/linear_elastic.cc b/deal.II/examples/step-41/source/linear_elastic.cc
deleted file mode 100644 (file)
index f90ad9b..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-#include "../include/linear_elastic.templates.h"
-
-
-template class LinearElastic<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/local_assemble.cc b/deal.II/examples/step-41/source/local_assemble.cc
deleted file mode 100644 (file)
index 3383328..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-#include "../include/local_assemble.templates.h"
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
-
-template class LocalAssemble<deal_II_dimension, DoFHandler<deal_II_dimension> >;
-template class LocalAssemble<deal_II_dimension, MGDoFHandler<deal_II_dimension> >;
diff --git a/deal.II/examples/step-41/source/local_assemble_base.cc b/deal.II/examples/step-41/source/local_assemble_base.cc
deleted file mode 100644 (file)
index 5ab1fc3..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-#include "../include/local_assemble_base.templates.h"
-
-template class LocalAssembleBase<deal_II_dimension, DoFHandler<deal_II_dimension> >;
-template class LocalAssembleBase<deal_II_dimension, MGDoFHandler<deal_II_dimension> >;
-
diff --git a/deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc b/deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc
deleted file mode 100644 (file)
index 1b704c5..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-#include "../include/local_assemble_elastic_matrix.templates.h"
-
-
-template class LocalAssembleElasticMatrix<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc b/deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc
deleted file mode 100644 (file)
index 0c351d4..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-#include "../include/local_assemble_elastic_rhs.templates.h"
-
-
-template class LocalAssembleElasticRHS<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/local_assemble_mass.cc b/deal.II/examples/step-41/source/local_assemble_mass.cc
deleted file mode 100644 (file)
index f6d6e9f..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-#include "../include/local_assemble_mass.templates.h"
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
-
-template class LocalAssembleMass<deal_II_dimension, DoFHandler<deal_II_dimension> >;
-template class LocalAssembleMass<deal_II_dimension, MGDoFHandler<deal_II_dimension> >;
diff --git a/deal.II/examples/step-41/source/local_assemble_plastic_project.cc b/deal.II/examples/step-41/source/local_assemble_plastic_project.cc
deleted file mode 100644 (file)
index a7600a7..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-#include "../include/local_assemble_plastic_project.templates.h"
-
-
-template class LocalAssemblePlasticProject<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/local_assemble_scalar_project.cc b/deal.II/examples/step-41/source/local_assemble_scalar_project.cc
deleted file mode 100644 (file)
index c8075e6..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-#include "../include/local_assemble_scalar_project.templates.h"
-
-
-template class LocalAssembleScalarProject<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/main.cc b/deal.II/examples/step-41/source/main.cc
deleted file mode 100644 (file)
index e45e1b5..0000000
+++ /dev/null
@@ -1,37 +0,0 @@
-#include <iostream>
-#include "../include/base.h"
-
-int main () 
-{
-  try
-      {
-        Base<deal_II_dimension> problem;
-        problem.run ();
-      }
-  catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-               << exc.what() << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      
-      return 1;
-    }
-  catch (...) 
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      return 1;
-    }
-
-  return 0;
-}
diff --git a/deal.II/examples/step-41/source/output_processor.cc b/deal.II/examples/step-41/source/output_processor.cc
deleted file mode 100644 (file)
index 2b28d33..0000000
+++ /dev/null
@@ -1,26 +0,0 @@
-#include "../include/output_processor.templates.h"
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/petsc_vector.h>
-#include <deal.II/lac/vector.h>
-
-template <>
-void OutputProcessor<deal_II_dimension,Vector<double> >::dump_vector (const Vector<double> &rhs, 
-                                                                   const std::string &filename )
-{
-    deallog.push("Dump");
-    deallog << "Writing: " << filename << std::endl;
-    std::ofstream out ( (filename).c_str());
-    rhs.block_write(out);
-    deallog.pop();
-}
-
-template class FilteredDataOut<deal_II_dimension, DoFHandler<deal_II_dimension> >;
-
-template class OutputProcessor<deal_II_dimension, Vector<double> >;
-template class OutputProcessor<deal_II_dimension, BlockVector<double> >;
-
-#ifdef DEAL_II_USE_PETSC
-template class OutputProcessor<deal_II_dimension, PETScWrappers::Vector >;
-#endif
-
-// template class OutputProcessor<deal_II_dimension, PETScWrappers::BlockVector >;
diff --git a/deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc b/deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc
deleted file mode 100644 (file)
index 3c1df05..0000000
+++ /dev/null
@@ -1,158 +0,0 @@
-/*
- * Immersed Boundary Problem:
- * 
- * Header Files
- *
- * Author: 
- * Luca Heltai <heltai@dimat.unipv.it>
- * =============
- * License: GPL.
- * =============
- * $Id: ibm_rhs.cc,v 1.34 2005/04/05 14:46:49 luca Exp $
- *
- */
-#include "../include/parsed_symmetric_tensor_function.templates.h"
-
-#if deal_II_dimension == 1
-template <>
-const SymmetricTensor<2, 1> & 
-  ParsedSymmetricTensorFunction<2, 1>::operator()(const Point<1> &p) const
-{ 
-  t[0][0] = f.value(p,0);
-  return t;
-}
-
-
-template <>
-const SymmetricTensor<4, 1> & 
-  ParsedSymmetricTensorFunction<4, 1>::operator()(const Point<1> &p) const
-{ 
-  t[0][0][0][0] = f.value(p, 0);
-  return t;
-}
-
-#endif
-
-#if deal_II_dimension == 2
-template <>
-const SymmetricTensor<2, 2> & 
-  ParsedSymmetricTensorFunction<2, 2>::operator()(const Point<2> &p) const
-{ 
-  t[0][0] = f.value(p,0);
-  t[0][1] = f.value(p,1);
-  t[1][1] = f.value(p,2);
-  return t;
-}
-
-
-template <>
-const SymmetricTensor<4, 2> & 
-  ParsedSymmetricTensorFunction<4, 2>::operator()(const Point<2> &p) const
-{ 
-  t[0][0][0][0] = f.value(p, 0);
-  
-  t[0][0][1][1] = f.value(p, 1);
-  t[1][1][0][0] = f.value(p, 1);
-  
-  t[0][0][0][1] = f.value(p, 2);
-  t[0][1][0][0] = f.value(p, 2);
-
-  t[1][1][1][1] = f.value(p, 3);
-  
-  t[1][1][0][1] = f.value(p, 4);
-  t[0][1][1][1] = f.value(p, 4);
-  
-  t[0][1][0][1] = f.value(p, 5);
-  t[1][0][0][1] = f.value(p, 5);
-  return t;
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-const SymmetricTensor<2, 3> & 
-  ParsedSymmetricTensorFunction<2, 3>::operator()(const Point<3> &p) const
-{ 
-  t[0][0] = f.value(p,0);
-  t[0][1] = f.value(p,1);
-  t[0][2] = f.value(p,2);
-  
-  t[1][1] = f.value(p,3);
-  t[1][2] = f.value(p,4);
-
-  t[2][2] = f.value(p,5);
-
-  return t;
-}
-
-
-template <>
-const SymmetricTensor<4, 3> & 
-  ParsedSymmetricTensorFunction<4, 3>::operator()(const Point<3> &p) const
-{ 
-  t[0][0][0][0] = f.value(p, 0);
-  
-  t[0][0][1][1] = f.value(p, 1);
-  t[1][1][0][0] = f.value(p, 1);
-  
-  t[0][0][2][2] = f.value(p, 2);
-  t[2][2][0][0] = f.value(p, 2);
-  
-  t[0][0][0][1] = f.value(p, 3);
-  t[0][1][0][0] = f.value(p, 3);
-  
-  t[0][0][1][2] = f.value(p, 4);
-  t[1][2][0][0] = f.value(p, 4);
-  
-  t[0][0][2][0] = f.value(p, 5);
-  t[2][0][0][0] = f.value(p, 5);
-  
-  t[1][1][1][1] = f.value(p, 6);
-  
-  t[1][1][2][2] = f.value(p, 7);
-  t[2][2][1][1] = f.value(p, 7);
-  
-  t[1][1][0][1] = f.value(p, 8);
-  t[0][1][1][1] = f.value(p, 8);
-  
-  t[1][1][1][2] = f.value(p, 9);
-  t[1][2][1][1] = f.value(p, 9);
-  
-  t[1][1][2][0] = f.value(p, 10);
-  t[2][0][1][1] = f.value(p, 10);
-  
-  t[2][2][2][2] = f.value(p, 11);
-  
-  t[2][2][0][1] = f.value(p, 12);
-  t[0][1][2][2] = f.value(p, 12);
-  
-  t[2][2][1][2] = f.value(p, 13);
-  t[1][2][2][2] = f.value(p, 13);
-  
-  t[2][2][0][2] = f.value(p, 14);
-  t[0][2][2][2] = f.value(p, 14);
-  
-  t[0][1][0][1] = f.value(p, 15);
-  
-  t[0][1][1][2] = f.value(p, 16);
-  t[1][2][0][1] = f.value(p, 16);
-  
-  t[0][1][0][2] = f.value(p, 17);
-  t[0][2][0][1] = f.value(p, 17);
-  
-  t[1][2][1][2] = f.value(p, 18);
-  
-  t[1][2][0][2] = f.value(p, 19);
-  t[0][2][1][2] = f.value(p, 19);
-  
-  t[0][2][0][2] = f.value(p, 20);
-
-  return t;
-}
-
-#endif
-
-template class ParsedSymmetricTensorFunction<2,deal_II_dimension>;
-template class ParsedSymmetricTensorFunction<4,deal_II_dimension>;
diff --git a/deal.II/examples/step-41/source/vector_space.cc b/deal.II/examples/step-41/source/vector_space.cc
deleted file mode 100644 (file)
index fa06236..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-#include "../include/vector_space.templates.h"
-
-template class VectorSpace<deal_II_dimension>;
diff --git a/deal.II/examples/step-41/step-41.cc b/deal.II/examples/step-41/step-41.cc
deleted file mode 100644 (file)
index e45e1b5..0000000
+++ /dev/null
@@ -1,37 +0,0 @@
-#include <iostream>
-#include "../include/base.h"
-
-int main () 
-{
-  try
-      {
-        Base<deal_II_dimension> problem;
-        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;
-}
index 657a12ade466d455070d5a060d38ede95175c7de..5401616f01a4b6c76332f7af15bb722165b77b63 100644 (file)
                                 // been covered in previous examples
                                 // and will thus not be further
                                 // commented on.
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/compressed_sparsity_pattern.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/numerics/data_out.h>
-
+#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/compressed_sparsity_pattern.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_gmres.h>
+#include <lac/precondition.h>
+#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 <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include "constrained_matrix.h"
 #include <fstream>
 #include <iostream>
 
@@ -51,7 +53,7 @@
                                 // tri-linear elements, but we will
                                 // now use it for bi-quadratic
                                 // elements:
-#include <deal.II/fe/fe_q.h>
+#include <fe/fe_q.h>
                                 // We will not read the grid from a
                                 // file as in the previous example,
                                 // but generate it using a function
@@ -61,7 +63,7 @@
                                 // the solution) in each step, so we
                                 // need the following include file
                                 // instead of <code>grid_in.h</code>:
-#include <deal.II/grid/grid_out.h>
+#include <grid/grid_out.h>
 
 
                                 // When using locally refined grids,
@@ -77,7 +79,7 @@
                                 // following file contains a class
                                 // which is used to handle these
                                 // constraints:
-#include <deal.II/lac/constraint_matrix.h>
+#include <lac/constraint_matrix.h>
 
                                 // In order to refine our grids
                                 // locally, we need a function from
@@ -86,7 +88,7 @@
                                 // coarsening based on the error
                                 // indicators we have computed. This
                                 // function is defined here:
-#include <deal.II/grid/grid_refinement.h>
+#include <grid/grid_refinement.h>
 
                                 // Finally, we need a simple way to
                                 // actually compute the refinement
                                 // often yields quite nicely adapted
                                 // grids for a wide class of
                                 // problems.
-#include <deal.II/numerics/error_estimator.h>
+#include <numerics/error_estimator.h>
 
                                 // Finally, this is as in previous
                                 // programs:
@@ -615,8 +617,9 @@ void Step6<dim>::assemble_system ()
                                   // constraints due to hanging
                                   // nodes. This is done using the
                                   // following two function calls:
-  hanging_node_constraints.condense (system_matrix);
-  hanging_node_constraints.condense (system_rhs);
+  // hanging_node_constraints.condense (system_matrix);
+  // hanging_node_constraints.condense (system_rhs);
+  
                                   // Using them, degrees of freedom
                                   // associated to hanging nodes have
                                   // been removed from the linear
@@ -695,16 +698,28 @@ void Step6<dim>::assemble_system ()
 template <int dim>
 void Step6<dim>::solve ()
 {
-  SolverControl           solver_control (1000, 1e-12);
-  SolverCG<>              solver (solver_control);
+  SolverControl           solver_control (10000, 1e-12);
+  SolverGMRES<>              solver (solver_control);
 
-  PreconditionSSOR<> preconditioner;
-  preconditioner.initialize(system_matrix, 1.2);
+  PreconditionIdentity preconditioner;
+//preconditioner.initialize(system_matrix, 1.2);rGM
 
-  solver.solve (system_matrix, solution, system_rhs,
-               preconditioner);
+  double scaling = 10000;
+  typename DoFHandler<dim>::active_cell_iterator
+    cell=dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for(; cell != endc; ++cell)
+    scaling = std::min(scaling, cell->diameter()*cell->diameter());
+
+  ConstrainedOperator<Vector<double>, SparseMatrix<double> >
+    cc(system_matrix, hanging_node_constraints, scaling);
 
-  hanging_node_constraints.distribute (solution);
+                                  // cc.distribute_rhs(system_rhs);
+  solution = system_rhs;
+  cc.distribute_rhs(system_rhs);
+  
+  solver.solve (cc, solution, system_rhs,
+               preconditioner);
 }
 
 
@@ -978,6 +993,19 @@ void Step6<dim>::output_results (const unsigned int cycle) const
 
   GridOut grid_out;
   grid_out.write_eps (triangulation, output);
+  
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  
+  filename = "solution-";
+  filename += ('0' + cycle);
+  filename += ".vtu";
+  output.close();
+  output.open(filename.c_str());
+  data_out.write_vtu (output);
 }
 
 
@@ -1090,8 +1118,8 @@ void Step6<dim>::run ()
   data_out.add_data_vector (solution, "solution");
   data_out.build_patches ();
 
-  std::ofstream output ("final-solution.eps");
-  data_out.write_eps (output);
+  std::ofstream output ("final-solution.vtk");
+  data_out.write_vtk (output);
 }
 
 

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.