]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Skeleton for quasistatic plasticity example: step-41
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 31 Aug 2010 14:24:46 +0000 (14:24 +0000)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 31 Aug 2010 14:24:46 +0000 (14:24 +0000)
git-svn-id: https://svn.dealii.org/trunk@21809 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/examples/step-41/CamClay.prm b/deal.II/examples/step-41/CamClay.prm
new file mode 100644 (file)
index 0000000..404fbd4
--- /dev/null
@@ -0,0 +1,306 @@
+# 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
new file mode 100644 (file)
index 0000000..7fb1643
--- /dev/null
@@ -0,0 +1,136 @@
+# $Id: Makefile.large 18849 2009-05-15 17:59:46Z bangerth $
+
+# 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 = 
+
+# 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
+
+
+# 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/base/include/base/*.h \
+                $D/lac/include/lac/*.h             \
+                $D/deal.II/include/*/*.h)
+
+# As before, define a list of libraries. This, of course depends on
+# the dimension in which we are working:
+libs.g   = $(lib-deal2-$(deal_II_dimension)d.g) \
+          $(lib-lac.g)                         \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-$(deal_II_dimension)d.o) \
+          $(lib-lac.o)                         \
+           $(lib-base.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
+
diff --git a/deal.II/examples/step-41/doc/builds-on b/deal.II/examples/step-41/doc/builds-on
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+
diff --git a/deal.II/examples/step-41/doc/intro.dox b/deal.II/examples/step-41/doc/intro.dox
new file mode 100644 (file)
index 0000000..446bf8b
--- /dev/null
@@ -0,0 +1,4 @@
+<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
new file mode 100644 (file)
index 0000000..56e049c
--- /dev/null
@@ -0,0 +1 @@
+solids
diff --git a/deal.II/examples/step-41/doc/results.dox b/deal.II/examples/step-41/doc/results.dox
new file mode 100644 (file)
index 0000000..6121fb6
--- /dev/null
@@ -0,0 +1,3 @@
+<h1>Results</h1>
+
+Astonishing results here.
diff --git a/deal.II/examples/step-41/doc/tooltip b/deal.II/examples/step-41/doc/tooltip
new file mode 100644 (file)
index 0000000..03e9264
--- /dev/null
@@ -0,0 +1 @@
+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
new file mode 100644 (file)
index 0000000..272ddb7
--- /dev/null
@@ -0,0 +1,122 @@
+//---------------------------------------------------------------------------
+//    $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 <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
new file mode 100644 (file)
index 0000000..05ad294
--- /dev/null
@@ -0,0 +1,113 @@
+#ifndef BASE_H
+#define BASE_H
+
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <base/parsed_function.h>
+#include <base/timer.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_generator.h>
+#include <numerics/data_out.h>
+
+#include <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
new file mode 100644 (file)
index 0000000..8664990
--- /dev/null
@@ -0,0 +1,468 @@
+#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.n_quadrature_points,
+                                                         std::vector<Tensor<1,dim> >(dim));
+  std::vector< Point<dim> > points(qf.n_quadrature_points);
+
+  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.n_quadrature_points; ++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.n_quadrature_points,
+                                                           std::vector<Tensor<1,dim> >(dim));
+  std::vector< std::vector< Tensor<1,dim> > > plastic_grads(qf.n_quadrature_points,
+                                                           std::vector<Tensor<1,dim> >(dim));
+  std::vector< Vector<double> > hard_values(qf.n_quadrature_points, Vector<double>(dim));
+  std::vector< Point<dim> > points(qf.n_quadrature_points);
+
+  //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
new file mode 100644 (file)
index 0000000..e926280
--- /dev/null
@@ -0,0 +1,172 @@
+#ifndef CAMCLAY
+#define CAMCLAY
+
+//deal.ii packages
+#include <base/logstream.h>
+#include <base/parameter_handler.h>
+#include <base/parsed_function.h>
+#include <base/symmetric_tensor.h>
+#include <base/table.h>
+#include <base/quadrature_lib.h>
+#include <base/utilities.h>
+#include <lac/sparse_direct.h>
+#include <lac/sparse_matrix.h>
+#include <lac/full_matrix.h>
+#include <lac/solver_control.h>
+#include <lac/solver_gmres.h>
+#include <fe/fe.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+#include <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
new file mode 100644 (file)
index 0000000..da3c493
--- /dev/null
@@ -0,0 +1,670 @@
+#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.n_quadrature_points;
+  
+
+  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.n_quadrature_points;
+
+  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
new file mode 100644 (file)
index 0000000..c8a2064
--- /dev/null
@@ -0,0 +1,72 @@
+#ifndef DOMAIN_H
+#define DOMAIN_H
+#include <fstream>
+
+#include <grid/grid_refinement.h>
+#include <grid/tria.h>
+#include <grid/grid_out.h>
+#include <base/parameter_handler.h>
+#include <base/path_search.h>
+#include <fe/mapping.h>
+#include <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
new file mode 100644 (file)
index 0000000..a1663ac
--- /dev/null
@@ -0,0 +1,120 @@
+#include "domain.h"
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_in.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_tools.h>
+#include <numerics/vectors.h>
+#include <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
new file mode 100644 (file)
index 0000000..239a251
--- /dev/null
@@ -0,0 +1,177 @@
+#ifndef LH_ERROR_HANDLER_H
+#define LH_ERROR_HANDLER_H
+
+#include <fstream>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <lac/vector.h>
+
+#include <grid/tria.h>
+
+// #include <numerics/error_estimator.h>
+#include <base/function.h>
+#include <numerics/solution_transfer.h>
+#include <numerics/data_out.h>
+
+#include <base/convergence_table.h>
+#include <base/logstream.h>
+#include <base/config.h>
+
+#include <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
new file mode 100644 (file)
index 0000000..93816d8
--- /dev/null
@@ -0,0 +1,396 @@
+#include "../include/error_handler.h"
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <base/utilities.h>
+
+#include <grid/grid_tools.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/mapping_q.h>
+#include <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
new file mode 100644 (file)
index 0000000..0d54b43
--- /dev/null
@@ -0,0 +1,77 @@
+#ifndef LINEARELASTIC
+#define LINEARELASTIC
+
+#include <base/logstream.h>
+#include <base/parameter_handler.h>
+#include <base/parsed_function.h>
+#include <lac/vector.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+#include <lac/filtered_matrix.h>
+#include <lac/precondition.h>
+#include <lac/sparse_ilu.h>
+#include <lac/solver_cg.h>
+#include <lac/sparse_direct.h>
+#include <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
new file mode 100644 (file)
index 0000000..02b51cf
--- /dev/null
@@ -0,0 +1,160 @@
+#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
new file mode 100644 (file)
index 0000000..26f5cf6
--- /dev/null
@@ -0,0 +1,66 @@
+#ifndef LOCAL_ASSEMBLE_STANDARD
+#define LOCAL_ASSEMBLE_STANDARD
+
+#include "local_assemble_base.h"
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <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
new file mode 100644 (file)
index 0000000..1a2d4a5
--- /dev/null
@@ -0,0 +1,67 @@
+#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
new file mode 100644 (file)
index 0000000..6e42d4d
--- /dev/null
@@ -0,0 +1,102 @@
+#ifndef __deal2__local_assemble_base_h
+#define __deal2__local_assemble_base_h
+
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <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
new file mode 100644 (file)
index 0000000..98fdc4e
--- /dev/null
@@ -0,0 +1,68 @@
+#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
new file mode 100644 (file)
index 0000000..41c096a
--- /dev/null
@@ -0,0 +1,50 @@
+#ifndef LOCAL_ASSEMBLE_ELASTIC_MATRIX
+#define LOCAL_ASSEMBLE_ELASTIC_MATRIX
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <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
new file mode 100644 (file)
index 0000000..3841fe3
--- /dev/null
@@ -0,0 +1,85 @@
+#include "../include/local_assemble_elastic_matrix.h"
+#include <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
new file mode 100644 (file)
index 0000000..d8bbbf3
--- /dev/null
@@ -0,0 +1,65 @@
+#ifndef LOCAL_ASSEMBLE_ELASTIC_RHS
+#define LOCAL_ASSEMBLE_ELASTIC_RHS
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <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
new file mode 100644 (file)
index 0000000..8954b90
--- /dev/null
@@ -0,0 +1,133 @@
+#include "../include/local_assemble_elastic_rhs.h"
+#include <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
new file mode 100644 (file)
index 0000000..4fa1de4
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef LOCAL_ASSEMBLE_HARDENING_PROJECT
+#define LOCAL_ASSEMBLE_HARDENING_PROJECT
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <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
new file mode 100644 (file)
index 0000000..501d967
--- /dev/null
@@ -0,0 +1,45 @@
+#ifndef LOCAL_ASSEMBLE_MASS_MATRIX
+#define LOCAL_ASSEMBLE_MASS_MATRIX
+
+#include "local_assemble.h"
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <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
new file mode 100644 (file)
index 0000000..cf33a56
--- /dev/null
@@ -0,0 +1,93 @@
+#include "../include/local_assemble_mass.h"
+#include "utilities.h"
+#include <base/quadrature_lib.h>
+#include <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
new file mode 100644 (file)
index 0000000..4a8710c
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef LOCAL_ASSEMBLE_PLASTIC_PROJECT
+#define LOCAL_ASSEMBLE_PLASTIC_PROJECT
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <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
new file mode 100644 (file)
index 0000000..4b74e2d
--- /dev/null
@@ -0,0 +1,93 @@
+#include "../include/local_assemble_plastic_project.h"
+#include <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
new file mode 100644 (file)
index 0000000..baf9515
--- /dev/null
@@ -0,0 +1,51 @@
+#ifndef LOCAL_ASSEMBLE_SCALAR_PROJECT
+#define LOCAL_ASSEMBLE_SCALAR_PROJECT
+
+#include "local_assemble.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <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
new file mode 100644 (file)
index 0000000..63118f3
--- /dev/null
@@ -0,0 +1,84 @@
+#include "../include/local_assemble_scalar_project.h"
+#include <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
new file mode 100644 (file)
index 0000000..2279632
--- /dev/null
@@ -0,0 +1,989 @@
+//---------------------------------------------------------------------------
+//    $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 <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <multigrid/mg_base.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <fe/fe.h>
+#include <base/quadrature_lib.h>
+
+#include <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
new file mode 100644 (file)
index 0000000..a37e80c
--- /dev/null
@@ -0,0 +1,176 @@
+#ifndef LH_OUTPUT_PROCESSOR_H
+#define LH_OUTPUT_PROCESSOR_H
+
+#include <fstream>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <lac/block_vector.h>
+
+#include <grid/grid_refinement.h>
+#include <grid/tria.h>
+#include <grid/filtered_iterator.h>
+
+// #include <numerics/error_estimator.h>
+#include <base/function.h>
+#include <numerics/solution_transfer.h>
+#include <numerics/data_out.h>
+
+#include <base/convergence_table.h>
+#include <base/logstream.h>
+#include <base/config.h>
+
+#include <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
new file mode 100644 (file)
index 0000000..0654966
--- /dev/null
@@ -0,0 +1,135 @@
+#include "../include/output_processor.h"
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <base/utilities.h>
+
+#include <grid/grid_tools.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/mapping_q.h>
+#include <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
new file mode 100644 (file)
index 0000000..8daeae4
--- /dev/null
@@ -0,0 +1,193 @@
+/*
+ * 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 <base/parameter_handler.h>
+#include <base/function.h>
+#include <base/function_parser.h>
+#include <base/point.h>
+#include <base/tensor.h>
+#include <base/symmetric_tensor.h>
+
+#include <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
new file mode 100644 (file)
index 0000000..15dc166
--- /dev/null
@@ -0,0 +1,135 @@
+/*
+ * 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 <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
new file mode 100644 (file)
index 0000000..7ff0f0f
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef UTILITIES_HLT
+#define UTILITIES_HLT
+
+#include <base/utilities.h>
+#include <grid/grid_tools.h>
+#include <lac/vector_memory.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <lac/block_vector.h>
+#include <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
new file mode 100644 (file)
index 0000000..17b50f9
--- /dev/null
@@ -0,0 +1,325 @@
+#ifndef VECTOR_SPACE_H
+#define VECTOR_SPACE_H
+
+#include <base/parameter_handler.h>
+#include <fe/fe_tools.h>
+#include <grid/persistent_tria.h>
+#include <multigrid/mg_dof_handler.h>
+#include <base/smartpointer.h>
+#include <fe/mapping.h>
+#include <dofs/dof_constraints.h>
+#include <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
new file mode 100644 (file)
index 0000000..3b62445
--- /dev/null
@@ -0,0 +1,583 @@
+#include "../include/vector_space.h"
+#include "../include/utilities.h"
+#include <base/logstream.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_in.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <fe/mapping_cartesian.h>
+#include <fe/mapping_q.h>
+#include <fe/fe.h>
+#include <base/utilities.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_tools.h>
+#include <numerics/solution_transfer.h>
+#include <numerics/vectors.h>
+#include <grid/grid_tools.h>
+#include <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
new file mode 100644 (file)
index 0000000..23f8127
--- /dev/null
@@ -0,0 +1,16 @@
+# 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
new file mode 100644 (file)
index 0000000..9da975a
--- /dev/null
@@ -0,0 +1,3 @@
+#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
new file mode 100644 (file)
index 0000000..506ae0c
--- /dev/null
@@ -0,0 +1,3 @@
+#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
new file mode 100644 (file)
index 0000000..fa3bc74
--- /dev/null
@@ -0,0 +1,3 @@
+#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
new file mode 100644 (file)
index 0000000..0403a7c
--- /dev/null
@@ -0,0 +1,11 @@
+#include "../include/error_handler.templates.h"
+#include <lac/block_vector.h>
+#include <lac/petsc_vector.h>
+#include <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
new file mode 100644 (file)
index 0000000..f90ad9b
--- /dev/null
@@ -0,0 +1,4 @@
+#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
new file mode 100644 (file)
index 0000000..54ac66c
--- /dev/null
@@ -0,0 +1,6 @@
+#include "../include/local_assemble.templates.h"
+#include <dofs/dof_handler.h>
+#include <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
new file mode 100644 (file)
index 0000000..5ab1fc3
--- /dev/null
@@ -0,0 +1,5 @@
+#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
new file mode 100644 (file)
index 0000000..1b704c5
--- /dev/null
@@ -0,0 +1,4 @@
+#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
new file mode 100644 (file)
index 0000000..0c351d4
--- /dev/null
@@ -0,0 +1,4 @@
+#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
new file mode 100644 (file)
index 0000000..78d0dc9
--- /dev/null
@@ -0,0 +1,6 @@
+#include "../include/local_assemble_mass.templates.h"
+#include <dofs/dof_handler.h>
+#include <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
new file mode 100644 (file)
index 0000000..a7600a7
--- /dev/null
@@ -0,0 +1,4 @@
+#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
new file mode 100644 (file)
index 0000000..c8075e6
--- /dev/null
@@ -0,0 +1,4 @@
+#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
new file mode 100644 (file)
index 0000000..e45e1b5
--- /dev/null
@@ -0,0 +1,37 @@
+#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
new file mode 100644 (file)
index 0000000..dd5ff00
--- /dev/null
@@ -0,0 +1,26 @@
+#include "../include/output_processor.templates.h"
+#include <lac/block_vector.h>
+#include <lac/petsc_vector.h>
+#include <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
new file mode 100644 (file)
index 0000000..3c1df05
--- /dev/null
@@ -0,0 +1,158 @@
+/*
+ * 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
new file mode 100644 (file)
index 0000000..fa06236
--- /dev/null
@@ -0,0 +1,3 @@
+#include "../include/vector_space.templates.h"
+
+template class VectorSpace<deal_II_dimension>;

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.