From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Tue, 20 May 2014 19:05:42 +0000 (+0000)
Subject: merge any_data branch
X-Git-Tag: v8.2.0-rc1~473
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c242b71a31a7cf8623a3152c01df7233261ab773;p=dealii.git

merge any_data branch

git-svn-id: https://svn.dealii.org/trunk@32947 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.h b/deal.II/include/deal.II/algorithms/theta_timestepping.h
index e6bb52410a..0fc7df05cf 100644
--- a/deal.II/include/deal.II/algorithms/theta_timestepping.h
+++ b/deal.II/include/deal.II/algorithms/theta_timestepping.h
@@ -114,7 +114,7 @@ namespace Algorithms
    * communicate the current time step information through AnyData as
    * well. Therefore, the AnyData objects handed as input to
    * #op_explicit and #op_implicit contain two entries of type
-   * <tt>const double</tt> reference named "Time" and "Timestep". Note
+   * <tt>const double*</tt> named "Time" and "Timestep". Note
    * that "Time" refers to the time at the beginning of the current
    * step for #op_explicit and at the end for #op_implicit,
    * respectively.
diff --git a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h
index 4072b37b4e..ab0d36a243 100644
--- a/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h
+++ b/deal.II/include/deal.II/algorithms/theta_timestepping.templates.h
@@ -91,8 +91,8 @@ namespace Algorithms
     // timestep
     AnyData src1;
     src1.add<const VECTOR*>(&solution, "Previous iterate");
-    src1.add<const double&>(d_explicit.time, "Time");
-    src1.add<const double&>(d_explicit.step, "Timestep");
+    src1.add<const double*>(&d_explicit.time, "Time");
+    src1.add<const double*>(&d_explicit.step, "Timestep");
     src1.merge(in);
 
     AnyData src2;
@@ -102,8 +102,8 @@ namespace Algorithms
     out1.add<VECTOR*>(aux, "Solution");
     // The data provided to the inner solver
     src2.add<const VECTOR*>(aux, "Previous time");
-    src2.add<const double&>(d_implicit.time, "Time");
-    src2.add<const double&>(d_explicit.step, "Timestep");
+    src2.add<const double*>(&d_implicit.time, "Time");
+    src2.add<const double*>(&d_explicit.step, "Timestep");
     src2.merge(in);
 
     if (output != 0)
diff --git a/deal.II/include/deal.II/meshworker/functional.h b/deal.II/include/deal.II/meshworker/functional.h
index a651ba1acc..0205951f21 100644
--- a/deal.II/include/deal.II/meshworker/functional.h
+++ b/deal.II/include/deal.II/meshworker/functional.h
@@ -19,6 +19,7 @@
 #define __deal2__mesh_worker_functional_h
 
 #include <deal.II/base/named_data.h>
+#include <deal.II/algorithms/any_data.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/base/mg_level_object.h>
 #include <deal.II/lac/block_vector.h>
@@ -107,49 +108,28 @@ namespace MeshWorker
     {
     public:
       /**
-       * The data type for
-       * communicating the cell and
-       * face vectors.
-       */
-      typedef NamedData<BlockVector<number>*> DataVectors;
-
-      /**
-       * The initialization
-       * function, specifying the
-       * @p results vectors and
-       * whether face data should
-       * be collected separately.
+       * The initialization function, specifying the @p results
+       * vectors and whether face data should be collected separately.
        *
-       * @p results should contain
-       * two block vectors named
-       * "cells" and "faces" (the
-       * latter only if
-       * @p separate_faces is
-       * true). In each of the two,
-       * each block should have
-       * equal size and be large
-       * enough to accommodate all
-       * user indices set in the
-       * cells and faces covered by
-       * the loop it is used
-       * in. Typically, for
-       * estimators, this is
-       * Triangulation::n_active_cells()
-       * and
-       * Triangulation::n_faces(),
-       * respectively.
+       * @p results should contain two block vectors named "cells" and
+       * "faces" (the latter only if @p separate_faces is true). In
+       * each of the two, each block should have equal size and be
+       * large enough to accommodate all user indices set in the cells
+       * and faces covered by the loop it is used in. Typically, for
+       * estimators, this is Triangulation::n_active_cells() and
+       * Triangulation::n_faces(), respectively.
        *
-       * The use of BlockVector may
-       * seem cumbersome, but it
-       * allows us to assemble
-       * several functionals at the
-       * same time, one in each
-       * block. The typical
-       * situation for error
-       * estimate is just having a
-       * single block in each vector.
+       * The use of BlockVector may seem cumbersome, but it allows us
+       * to assemble several functionals at the same time, one in each
+       * block. The typical situation for error estimate is just
+       * having a single block in each vector.
        */
-      void initialize(DataVectors &results,
+      void initialize(AnyData &results, bool separate_faces = true);
+
+	/**
+	 * @deprecated
+	 */
+	void initialize(NamedData<BlockVector<number>*> &results,
                       bool separate_faces = true);
       /**
        * Initialize the local data
@@ -188,7 +168,7 @@ namespace MeshWorker
        */
       number operator() (const unsigned int i) const;
     private:
-      DataVectors results;
+      AnyData results;
       bool separate_faces;
     };
 //----------------------------------------------------------------------//
@@ -247,13 +227,29 @@ namespace MeshWorker
 
     template <typename number>
     inline void
-    CellsAndFaces<number>::initialize(DataVectors &r, bool sep)
+    CellsAndFaces<number>::initialize(AnyData &r, bool sep)
+    {
+      Assert(r.name(0) == "cells", AnyData::ExcNameMismatch(0, "cells"));
+      if (sep)
+        {
+          Assert(r.name(1) == "faces", AnyData::ExcNameMismatch(1, "faces"));
+          AssertDimension(r.entry<BlockVector<double>*>(0)->n_blocks(),
+			  r.entry<BlockVector<double>*>(1)->n_blocks());
+        }
+
+      results = r;
+      separate_faces = sep;
+    }
+
+    template <typename number>
+    inline void
+    CellsAndFaces<number>::initialize(NamedData<BlockVector<number>*> &r, bool sep)
     {
-      Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells"));
+      Assert(r.name(0) == "cells", AnyData::ExcNameMismatch(0, "cells"));
       if (sep)
         {
-          Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces"));
-          AssertDimension(r(0)->n_blocks(), r(1)->n_blocks());
+          Assert(r.name(1) == "faces", AnyData::ExcNameMismatch(1, "faces"));
+          AssertDimension(r(0)->n_blocks(),r(1)->n_blocks());
         }
 
       results = r;
@@ -266,7 +262,7 @@ namespace MeshWorker
     inline void
     CellsAndFaces<number>::initialize_info(DOFINFO &info, bool) const
     {
-      info.initialize_numbers(results(0)->n_blocks());
+      info.initialize_numbers(results.entry<BlockVector<double>*>(0)->n_blocks());
     }
 
 
@@ -275,14 +271,15 @@ namespace MeshWorker
     inline void
     CellsAndFaces<number>::assemble(const DOFINFO &info)
     {
+      BlockVector<double>* v;
+      if (separate_faces &&
+	  info.face_number != deal_II_numbers::invalid_unsigned_int)
+	v = results.entry<BlockVector<double>*>(1);
+      else
+	v = results.entry<BlockVector<double>*>(0);
+      
       for (unsigned int i=0; i<info.n_values(); ++i)
-        {
-          if (separate_faces &&
-              info.face_number != deal_II_numbers::invalid_unsigned_int)
-            results(1)->block(i)(info.face->user_index()) += info.value(i);
-          else
-            results(0)->block(i)(info.cell->user_index()) += info.value(i);
-        }
+	v->block(i)(info.cell->user_index()) += info.value(i);
     }
 
 
@@ -296,15 +293,17 @@ namespace MeshWorker
         {
           if (separate_faces)
             {
+	      BlockVector<double>* v1 = results.entry<BlockVector<double>*>(1);
               const double J = info1.value(i) + info2.value(i);
-              results(1)->block(i)(info1.face->user_index()) += J;
+              v1->block(i)(info1.face->user_index()) += J;
               if (info2.face != info1.face)
-                results(1)->block(i)(info2.face->user_index()) += J;
+                v1->block(i)(info2.face->user_index()) += J;
             }
           else
             {
-              results(0)->block(i)(info1.cell->user_index()) += .5*info1.value(i);
-              results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i);
+	      BlockVector<double>* v0 = results.entry<BlockVector<double>*>(0);
+              v0->block(i)(info1.cell->user_index()) += .5*info1.value(i);
+              v0->block(i)(info2.cell->user_index()) += .5*info2.value(i);
             }
         }
     }
diff --git a/deal.II/include/deal.II/meshworker/integration_info.h b/deal.II/include/deal.II/meshworker/integration_info.h
index 4eb299711e..a120eced0f 100644
--- a/deal.II/include/deal.II/meshworker/integration_info.h
+++ b/deal.II/include/deal.II/meshworker/integration_info.h
@@ -248,7 +248,7 @@ namespace MeshWorker
    * undertaken to use this class.
    *
    * First, you should consider if you need values from any vectors in a
-   * NamedData object. If so, fill the VectorSelector objects
+   * AnyData object. If so, fill the VectorSelector objects
    * #cell_selector, #boundary_selector and #face_selector with their names
    * and the data type (value, gradient, Hessian) to be extracted.
    *
@@ -313,17 +313,22 @@ namespace MeshWorker
      * formulas, which integrate polynomial bilinear forms exactly.
      */
     template <typename VECTOR>
+    void initialize(const FiniteElement<dim, spacedim> &el,
+                    const Mapping<dim, spacedim> &mapping,
+                    const AnyData &data,
+		    const VECTOR& dummy,
+                    const BlockInfo *block_info = 0);
+    /**
+     * @deprecated Use AnyData instead of NamedData.
+     */
+    template <typename VECTOR>
     void initialize(const FiniteElement<dim, spacedim> &el,
                     const Mapping<dim, spacedim> &mapping,
                     const NamedData<VECTOR *> &data,
                     const BlockInfo *block_info = 0);
 
     /**
-     * Initialize the IntegrationInfo objects contained.
-     *
-     * Before doing so, add update flags necessary to produce the data
-     * needed and also set uninitialized quadrature rules to Gauss
-     * formulas, which integrate polynomial bilinear forms exactly.
+     * @deprecated Use AnyData instead of NamedData.
      */
     template <typename VECTOR>
     void initialize(const FiniteElement<dim, spacedim> &el,
@@ -767,6 +772,44 @@ namespace MeshWorker
   }
 
 
+  template <int dim, int sdim>
+  template <typename VECTOR>
+  void
+  IntegrationInfoBox<dim,sdim>::initialize(
+    const FiniteElement<dim,sdim> &el,
+    const Mapping<dim,sdim> &mapping,
+    const AnyData &data,
+    const VECTOR& dummy,
+    const BlockInfo *block_info)
+  {
+    initialize(el, mapping, block_info);
+    std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> > p;
+    VectorDataBase<dim,sdim>* pp;
+    
+    p = std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (cell_selector));
+    // Public member function of parent class was not found without
+    // explicit cast
+    pp = &*p;
+    pp->initialize(data);
+    cell_data = p;
+    cell.initialize_data(p);
+
+    p = std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (boundary_selector));
+    pp = &*p;
+    pp->initialize(data);
+    boundary_data = p;
+    boundary.initialize_data(p);
+
+    p = std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (face_selector));
+    pp = &*p;
+    pp->initialize(data);
+    face_data = p;
+    face.initialize_data(p);
+    subface.initialize_data(p);
+    neighbor.initialize_data(p);
+  }
+
+
   template <int dim, int sdim>
   template <typename VECTOR>
   void
diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h
index ed692c57f3..297e2421a6 100644
--- a/deal.II/include/deal.II/meshworker/simple.h
+++ b/deal.II/include/deal.II/meshworker/simple.h
@@ -19,6 +19,7 @@
 #define __deal2__mesh_worker_simple_h
 
 #include <deal.II/base/named_data.h>
+#include <deal.II/algorithms/any_data.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/base/mg_level_object.h>
 #include <deal.II/lac/block_vector.h>
@@ -60,8 +61,18 @@ namespace MeshWorker
     class ResidualSimple
     {
     public:
-      void initialize(NamedData<VECTOR *> &results);
+     /**
+      * Initialize with an AnyData object holding the result of
+      * assembling.
+      *
+      * Assembling currently writes into the first vector of <tt>results</tt>.
+      */
+      void initialize(AnyData& results);
 
+      /**
+       * @deprecated Use initialize(AnyData&) instead.
+       */
+      void initialize(NamedData<VECTOR *> &results);
       /**
        * Initialize the constraints.
        */
@@ -114,7 +125,7 @@ namespace MeshWorker
       /**
        * The global residal vectors filled by assemble().
        */
-      NamedData<SmartPointer<VECTOR,ResidualSimple<VECTOR> > > residuals;
+      AnyData residuals;
       /**
        * A pointer to the object containing constraints.
        */
@@ -500,6 +511,13 @@ namespace MeshWorker
 
 //----------------------------------------------------------------------//
 
+    template <class VECTOR>
+    inline void
+    ResidualSimple<VECTOR>::initialize(AnyData &results)
+    {
+      residuals = results;
+    }
+
     template <class VECTOR>
     inline void
     ResidualSimple<VECTOR>::initialize(NamedData<VECTOR *> &results)
@@ -537,18 +555,19 @@ namespace MeshWorker
     {
       for (unsigned int k=0; k<residuals.size(); ++k)
         {
+	  VECTOR* v = residuals.entry<VECTOR*>(k);
           if (constraints == 0)
             {
               for (unsigned int i=0; i<info.vector(k).block(0).size(); ++i)
-                (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
+                (*v)(info.indices[i]) += info.vector(k).block(0)(i);
             }
           else
             {
               if (info.indices_by_block.size() == 0)
-                constraints->distribute_local_to_global(info.vector(k).block(0), info.indices, (*residuals(k)));
+                constraints->distribute_local_to_global(info.vector(k).block(0), info.indices, *v);
               else
                 for (unsigned int i=0; i != info.vector(k).n_blocks(); ++i)
-                  constraints->distribute_local_to_global(info.vector(k).block(i), info.indices_by_block[i], (*residuals(k)));
+                  constraints->distribute_local_to_global(info.vector(k).block(i), info.indices_by_block[i], *v);
             }
         }
     }
@@ -561,30 +580,31 @@ namespace MeshWorker
     {
       for (unsigned int k=0; k<residuals.size(); ++k)
         {
+	  VECTOR* v = residuals.entry<VECTOR*>(k);
           if (constraints == 0)
             {
               for (unsigned int i=0; i<info1.vector(k).block(0).size(); ++i)
-                (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
+                (*v)(info1.indices[i]) += info1.vector(k).block(0)(i);
               for (unsigned int i=0; i<info2.vector(k).block(0).size(); ++i)
-                (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
+                (*v)(info2.indices[i]) += info2.vector(k).block(0)(i);
             }
           else
             {
               if (info1.indices_by_block.size() == 0 && info2.indices_by_block.size() == 0)
                 {
                   constraints->distribute_local_to_global
-                  (info1.vector(k).block(0), info1.indices, (*residuals(k)));
+                  (info1.vector(k).block(0), info1.indices, *v);
                   constraints->distribute_local_to_global
-                  (info2.vector(k).block(0), info2.indices, (*residuals(k)));
+                  (info2.vector(k).block(0), info2.indices, *v);
                 }
               else if (info1.indices_by_block.size() != 0 && info2.indices_by_block.size() != 0)
                 {
                   for (unsigned int i=0; i<info1.vector(k).n_blocks(); ++i)
                     {
                       constraints->distribute_local_to_global
-                      (info1.vector(k).block(i), info1.indices_by_block[i], (*residuals(k)));
+                      (info1.vector(k).block(i), info1.indices_by_block[i], *v);
                       constraints->distribute_local_to_global
-                      (info2.vector(k).block(i), info2.indices_by_block[i], (*residuals(k)));
+                      (info2.vector(k).block(i), info2.indices_by_block[i], *v);
                     }
                 }
               else
@@ -1132,7 +1152,7 @@ namespace MeshWorker
     inline void
     SystemSimple<MATRIX,VECTOR>::initialize(MATRIX &m, VECTOR &rhs)
     {
-      NamedData<VECTOR *> data;
+      AnyData data;
       VECTOR *p = &rhs;
       data.add(p, "right hand side");
 
diff --git a/deal.II/include/deal.II/meshworker/vector_selector.h b/deal.II/include/deal.II/meshworker/vector_selector.h
index 5976df4dbd..ff4d2e2fd4 100644
--- a/deal.II/include/deal.II/meshworker/vector_selector.h
+++ b/deal.II/include/deal.II/meshworker/vector_selector.h
@@ -18,6 +18,7 @@
 #define __deal2__mesh_worker_vector_selector_h
 
 #include <deal.II/base/named_data.h>
+#include <deal.II/algorithms/any_data.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/base/mg_level_object.h>
@@ -32,7 +33,7 @@ namespace MeshWorker
   /**
    * A class that selects vectors from a list of named vectors.
    *
-   * Since the number of vectors in a NamedData object may grow with every
+   * Since the number of vectors in an AnyData object may grow with every
    * nesting of applications or loops, it is important to be able to
    * select those, which are actually used in computing residuals etc.
    * This class organizes the selection.
@@ -49,23 +50,13 @@ namespace MeshWorker
   {
   public:
     /**
-     * Add a vector to the
-     * selection of finite element
-     * functions. The arguments are
-     * the name of the vector and
-     * indicators, which
-     * information is to be
-     * extracted from the
-     * vector. The name refers to
-     * an entry in a NamedData
-     * object, which will be
-     * identified by initialize().
-     * The three bool parameters
-     * indicate, whether values,
-     * greadients and Hessians of
-     * the finite element function
-     * are to be computed on each
-     * cell or face.
+     * Add a vector to the selection of finite element functions. The
+     * arguments are the name of the vector and indicators, which
+     * information is to be extracted from the vector. The name refers
+     * to an entry in a AnyData object, which will be identified by
+     * initialize().  The three bool parameters indicate, whether
+     * values, greadients and Hessians of the finite element function
+     * are to be computed on each cell or face.
      */
     void add(const std::string &name,
              const bool values = true,
@@ -73,10 +64,8 @@ namespace MeshWorker
              const bool hessians = false);
 
     /**
-     * Does the same as the function
-     * above but it is possible to
-     * select a block of the global
-     * vector.
+     * Does the same as the function above but it is possible to
+     * select a block of the global vector.
      */
 //      void add(const std::string& name,
 //               const unsigned int selected_block,
@@ -85,21 +74,19 @@ namespace MeshWorker
 //             bool hessians = false);
 
     /**
-     * Initialize the selection
-     * field with a data
-     * vector. While add() only
-     * enters the names of vectors,
-     * which will be used in the integration
-     * loop over cells and faces,
-     * this function links the
-     * names to actual vectos in a
-     * NamedData object.
+     * Initialize the selection field with a data vector. While add()
+     * only enters the names of vectors, which will be used in the
+     * integration loop over cells and faces, this function links the
+     * names to actual vectos in a AnyData object.
      *
-     * @note This function caches
-     * the index associated with a
-     * name. Therefore, it must be
-     * called every time after the
-     * NamedData object has changed.
+     * @note This function caches the index associated with a
+     * name. Therefore, it must be called every time after the AnyData
+     * object has changed.
+     */
+    void initialize(const AnyData &);
+
+    /**
+     * @deprecated Use AnyData instead of NamedData.
      */
     template <class DATA>
     void initialize(const NamedData<DATA> &);
@@ -110,20 +97,17 @@ namespace MeshWorker
     bool empty () const;
 
     /**
-     * Returns true if values are
-     * selected for any vector.
+     * Returns true if values are selected for any vector.
      */
     bool has_values () const;
 
     /**
-     * Returns true if gradients are
-     * selected for any vector.
+     * Returns true if gradients are selected for any vector.
      */
     bool has_gradients () const;
 
     /**
-     * Returns true if hessians are
-     * selected for any vector.
+     * Returns true if hessians are selected for any vector.
      */
     bool has_hessians () const;
 
@@ -162,6 +146,12 @@ namespace MeshWorker
      * selection to the stream.
      */
     template <class STREAM, typename DATA>
+    void print (STREAM &s, const AnyData &v) const;
+
+    /**
+     * @deprecated Use AnyData instead of NamedData.
+     */
+    template <class STREAM, typename DATA>
     void print (STREAM &s, const NamedData<DATA> &v) const;
 
     /**
@@ -183,14 +173,12 @@ namespace MeshWorker
     NamedSelection value_selection;
 
     /**
-     * Selection of the vectors
-     * used to compute gradients.
+     * Selection of the vectors used to compute gradients.
      */
     NamedSelection gradient_selection;
 
     /**
-     * Selection of the vectors
-     * used to compute hessians.
+     * Selection of the vectors used to compute hessians.
      */
     NamedSelection hessian_selection;
   };
@@ -213,81 +201,61 @@ namespace MeshWorker
     VectorDataBase();
 
     /**
-     * Constructor from a base
-     * class object
+     * Constructor from a base class object
      */
     VectorDataBase(const VectorSelector &);
 
     /**
-     * Virtual, but empty
-     * destructor.
+     * Initialize with a AnyData object and cache the indices in the
+     * VectorSelector base class.
+     *
+     * @note Make sure the VectorSelector base class was filled with
+     * reasonable data before calling this function.
+     */
+    void initialize(const AnyData&);
+
+    /**
+     * Virtual, but empty destructor.
      */
     virtual ~VectorDataBase();
     /**
-     * The only function added to
-     * VectorSelector is an
-     * abstract virtual function
-     * implemented in the derived
-     * class template and called by
-     * IntegrationInfo.
+     * The only function added to VectorSelector is an abstract
+     * virtual function implemented in the derived class template and
+     * called by IntegrationInfo.
      *
-     * Depending on the selections
-     * made in our base class, this
-     * fills the first three
-     * arguments with the local
-     * data of the finite element
-     * functions. It is usually
-     * called either for the whole
-     * FESystem, or for each base
-     * element separately.
+     * Depending on the selections made in our base class, this fills
+     * the first three arguments with the local data of the finite
+     * element functions. It is usually called either for the whole
+     * FESystem, or for each base element separately.
      *
-     * @param values is the vector
-     * filled with the values of
-     * the finite element function
-     * in the quadrature points.
+     * @param values is the vector filled with the values of the
+     * finite element function in the quadrature points.
      *
-     * @param gradients is the vector
-     * filled with the derivatives of
-     * the finite element function
-     * in the quadrature points.
+     * @param gradients is the vector filled with the derivatives of
+     * the finite element function in the quadrature points.
      *
-     * @param hessians is the
-     * vector filled with the
-     * second derivatives of the
-     * finite element function in
-     * the quadrature points.
+     * @param hessians is the vector filled with the second
+     * derivatives of the finite element function in the quadrature
+     * points.
      *
-     * @param fe is the
-     * FEValuesBase object which is
-     * used to compute the function
-     * values. Its UpdateFlags must
-     * have been set appropriately.
+     * @param fe is the FEValuesBase object which is used to compute
+     * the function values. Its UpdateFlags must have been set
+     * appropriately.
      *
-     * @param index is the local
-     * index vector. If @p fe
-     * refers to base elements of
-     * the system, this vector
-     * should be sorted by block
-     * and the arguments @p start
-     * and @p size below specify
-     * the subset of @p indices
-     * used.
+     * @param index is the local index vector. If @p fe refers to base
+     * elements of the system, this vector should be sorted by block
+     * and the arguments @p start and @p size below specify the subset
+     * of @p indices used.
      *
-     * @param component is the first
-     * index in @p values, @p
-     * gradients and @p hessians
-     * entered in this function.
+     * @param component is the first index in @p values, @p gradients
+     * and @p hessians entered in this function.
      *
-     * @param n_comp is the number
-     * of components to be filled.
+     * @param n_comp is the number of components to be filled.
      *
-     * @param start is the first
-     * index of this block in @p
-     * indices, or zero if
-     * no base elements are used.
+     * @param start is the first index of this block in @p indices, or
+     * zero if no base elements are used.
      *
-     * @param size is the number of
-     * dofs per cell of the current
+     * @param size is the number of dofs per cell of the current
      * element or base element.
      */
     virtual void fill(
@@ -302,15 +270,10 @@ namespace MeshWorker
       const unsigned int size) const;
 
     /**
-     * Fill the local data vector
-     * from level vectors. Performs
-     * exactly what the other
-     * fill() does, but uses the
-     * cell level to access a
-     * single level out of a
-     * hierarchy of level vectors,
-     * instead of a global data
-     * vector on the active cells.
+     * Fill the local data vector from level vectors. Performs exactly
+     * what the other fill() does, but uses the cell level to access a
+     * single level out of a hierarchy of level vectors, instead of a
+     * global data vector on the active cells.
      */
     virtual void mg_fill(
       std::vector<std::vector<std::vector<double> > > &values,
@@ -328,13 +291,15 @@ namespace MeshWorker
      * The memory used by this object.
      */
     std::size_t memory_consumption () const;
+  protected:
+    AnyData data;
   };
 
 
   /**
    * Based on VectorSelector, this is the class that implements the
    * function VectorDataBase::fill() for a certain type of vector, using
-   * NamedData to identify vectors by name.
+   * AnyData to identify vectors by name.
    *
    * @ingroup MeshWorker
    * @author Guido Kanschat, 2009
@@ -355,30 +320,16 @@ namespace MeshWorker
     VectorData(const VectorSelector &);
 
     /**
-     * Initialize with a NamedData
-     * object and cache the indices
-     * in the VectorSelector base
-     * class.
-     *
-     * @note Make sure the
-     * VectorSelector base class
-     * was filled with reasonable
-     * data before calling this
-     * function.
+     * @deprecated Use AnyData instead of NamedData.
      */
     void initialize(const NamedData<VECTOR *> &);
 
     /**
-     * Initialize with a single
-     * vector and cache the indices
-     * in the VectorSelector base
-     * class.
+     * Initialize with a single vector and cache the indices in the
+     * VectorSelector base class.
      *
-     * @note Make sure the
-     * VectorSelector base class
-     * was filled with reasonable
-     * data before calling this
-     * function.
+     * @note Make sure the VectorSelector base class was filled with
+     * reasonable data before calling this function.
      */
     void initialize(const VECTOR *, const std::string &name);
 
@@ -393,20 +344,29 @@ namespace MeshWorker
       const unsigned int start,
       const unsigned int size) const;
 
-    /**
+    virtual void mg_fill(
+      std::vector<std::vector<std::vector<double> > > &values,
+      std::vector<std::vector<std::vector<Tensor<1,dim> > > > &gradients,
+      std::vector<std::vector<std::vector<Tensor<2,dim> > > > &hessians,
+      const FEValuesBase<dim,spacedim> &fe,
+      const unsigned int level,
+      const std::vector<types::global_dof_index> &index,
+      const unsigned int component,
+      const unsigned int n_comp,
+      const unsigned int start,
+      const unsigned int size) const;
+
+     /**
      * The memory used by this object.
      */
     std::size_t memory_consumption () const;
-
-  private:
-    NamedData<SmartPointer<const VECTOR,VectorData<VECTOR,dim,spacedim> > > data;
   };
 
 
   /**
    * Based on VectorSelector, this is the class that implements the
    * function VectorDataBase::fill() for a certain type of multilevel vectors, using
-   * NamedData to identify vectors by name.
+   * AnyData to identify vectors by name.
    *
    * @ingroup MeshWorker
    * @author Guido Kanschat, 2010
@@ -421,58 +381,29 @@ namespace MeshWorker
      */
     MGVectorData();
     /**
-     * Constructor using a
-     * prefilled VectorSelector
+     * Constructor using a prefilled VectorSelector
      */
     MGVectorData(const VectorSelector &);
 
     /**
-     * Initialize with a NamedData
-     * object and cache the indices
-     * in the VectorSelector base
-     * class.
-     *
-     * @note Make sure the
-     * VectorSelector base class
-     * was filled with reasonable
-     * data before calling this
-     * function.
+     * @deprecated Use AnyData instead of NamedData.
      */
     void initialize(const NamedData<MGLevelObject<VECTOR>*> &);
 
     /**
-     * Initialize with a single
-     * vector and cache the indices
-     * in the VectorSelector base
-     * class.
+     * Initialize with a single vector and cache the indices in the
+     * VectorSelector base class.
      *
-     * @note Make sure the
-     * VectorSelector base class
-     * was filled with reasonable
-     * data before calling this
-     * function.
+     * @note Make sure the VectorSelector base class was filled with
+     * reasonable data before calling this function.
      */
     void initialize(const MGLevelObject<VECTOR> *, const std::string &name);
 
-    virtual void mg_fill(
-      std::vector<std::vector<std::vector<double> > > &values,
-      std::vector<std::vector<std::vector<Tensor<1,dim> > > > &gradients,
-      std::vector<std::vector<std::vector<Tensor<2,dim> > > > &hessians,
-      const FEValuesBase<dim,spacedim> &fe,
-      const unsigned int level,
-      const std::vector<types::global_dof_index> &index,
-      const unsigned int component,
-      const unsigned int n_comp,
-      const unsigned int start,
-      const unsigned int size) const;
 
     /**
      * The memory used by this object.
      */
     std::size_t memory_consumption () const;
-
-  private:
-    NamedData<SmartPointer<const MGLevelObject<VECTOR>,MGVectorData<VECTOR,dim,spacedim> > > data;
   };
 
 
@@ -504,6 +435,14 @@ namespace MeshWorker
   //}
 
 
+  inline void
+  VectorSelector::initialize(const AnyData &src)
+  {
+    value_selection.initialize(src);
+    gradient_selection.initialize(src);
+    hessian_selection.initialize(src);
+  }
+
   template <typename DATA>
   inline void
   VectorSelector::initialize(const NamedData<DATA> &src)
@@ -596,6 +535,23 @@ namespace MeshWorker
   }
 
 
+  template <class STREAM, typename DATA>
+  inline void
+  VectorSelector::print(STREAM &s, const AnyData &v) const
+  {
+    s << "values:   ";
+    for (unsigned int i=0; i<n_values(); ++i)
+      s << " '" << v.name(value_selection(i)) << '\'';
+    s << std::endl << "gradients:";
+    for (unsigned int i=0; i<n_gradients(); ++i)
+      s << " '" << v.name(gradient_selection(i)) << '\'';
+    s << std::endl << "hessians: ";
+    for (unsigned int i=0; i<n_hessians(); ++i)
+      s << " '" << v.name(hessian_selection(i)) << '\'';
+    s << std::endl;
+  }
+
+
   template <class STREAM, typename DATA>
   inline void
   VectorSelector::print(STREAM &s, const NamedData<DATA> &v) const
diff --git a/deal.II/include/deal.II/meshworker/vector_selector.templates.h b/deal.II/include/deal.II/meshworker/vector_selector.templates.h
index dc0e6fbbff..c48e19e25d 100644
--- a/deal.II/include/deal.II/meshworker/vector_selector.templates.h
+++ b/deal.II/include/deal.II/meshworker/vector_selector.templates.h
@@ -40,6 +40,15 @@ namespace MeshWorker
   {}
 
 
+  template <int dim, int spacedim>
+  void
+  VectorDataBase<dim, spacedim>::initialize(const AnyData &d)
+  {
+    this->data = d;
+    VectorSelector::initialize(d);
+  }
+
+
   template <int dim, int spacedim>
   void
   VectorDataBase<dim, spacedim>::fill(
@@ -93,7 +102,7 @@ namespace MeshWorker
   void
   VectorData<VECTOR, dim, spacedim>::initialize(const NamedData<VECTOR *> &d)
   {
-    data = d;
+    this->data = d;
     VectorSelector::initialize(d);
   }
 
@@ -103,8 +112,8 @@ namespace MeshWorker
   VectorData<VECTOR, dim, spacedim>::initialize(const VECTOR *v, const std::string &name)
   {
     SmartPointer<const VECTOR,VectorData<VECTOR, dim, spacedim> > p = v;
-    data.add(p, name);
-    VectorSelector::initialize(data);
+    this->data.add(p, name);
+    VectorSelector::initialize(this->data);
   }
 
 
@@ -125,25 +134,26 @@ namespace MeshWorker
     AssertDimension(gradients.size(), this->n_gradients());
     AssertDimension(hessians.size(), this->n_hessians());
 
+    const AnyData& data = this->data;
     for (unsigned int i=0; i<this->n_values(); ++i)
       {
-        const VECTOR &src = *data(this->value_index(i));
+        const VECTOR* src = data.read_ptr<VECTOR>(this->value_index(i));
         VectorSlice<std::vector<std::vector<double> > > dst(values[i], component, n_comp);
-        fe.get_function_values(src, make_slice(index, start, size), dst, true);
+        fe.get_function_values(*src, make_slice(index, start, size), dst, true);
       }
 
     for (unsigned int i=0; i<this->n_gradients(); ++i)
       {
-        const VECTOR &src = *data(this->gradient_index(i));
+        const VECTOR* src = data.read_ptr<VECTOR>(this->gradient_index(i));
         VectorSlice<std::vector<std::vector<Tensor<1,dim> > > > dst(gradients[i], component, n_comp);
-        fe.get_function_gradients(src, make_slice(index, start, size), dst, true);
+        fe.get_function_gradients(*src, make_slice(index, start, size), dst, true);
       }
 
     for (unsigned int i=0; i<this->n_hessians(); ++i)
       {
-        const VECTOR &src = *data(this->hessian_index(i));
+        const VECTOR* src = data.read_ptr<VECTOR>(this->hessian_index(i));
         VectorSlice<std::vector<std::vector<Tensor<2,dim> > > > dst(hessians[i], component, n_comp);
-        fe.get_function_hessians(src, make_slice(index, start, size), dst, true);
+        fe.get_function_hessians(*src, make_slice(index, start, size), dst, true);
       }
   }
 
@@ -153,7 +163,7 @@ namespace MeshWorker
   VectorData<VECTOR, dim, spacedim>::memory_consumption () const
   {
     std::size_t mem = VectorSelector::memory_consumption();
-    mem += sizeof (data);
+    mem += sizeof (this->data);
     return mem;
   }
 
@@ -175,7 +185,7 @@ namespace MeshWorker
   void
   MGVectorData<VECTOR, dim, spacedim>::initialize(const NamedData<MGLevelObject<VECTOR>*> &d)
   {
-    data = d;
+    this->data = d;
     VectorSelector::initialize(d);
   }
 
@@ -186,14 +196,14 @@ namespace MeshWorker
   {
     SmartPointer<const MGLevelObject<VECTOR>, MGVectorData<VECTOR, dim, spacedim> >
     p = v;
-    data.add(p, name);
-    VectorSelector::initialize(data);
+    this->data.add(p, name);
+    VectorSelector::initialize(this->data);
   }
 
 
   template <class VECTOR, int dim, int spacedim>
   void
-  MGVectorData<VECTOR, dim, spacedim>::mg_fill(
+  VectorData<VECTOR, dim, spacedim>::mg_fill(
     std::vector<std::vector<std::vector<double> > > &values,
     std::vector<std::vector<std::vector<Tensor<1,dim> > > > &gradients,
     std::vector<std::vector<std::vector<Tensor<2,dim> > > > &hessians,
@@ -209,25 +219,26 @@ namespace MeshWorker
     AssertDimension(gradients.size(), this->n_gradients());
     AssertDimension(hessians.size(), this->n_hessians());
 
+    const AnyData& data = this->data;
     for (unsigned int i=0; i<this->n_values(); ++i)
       {
-        const VECTOR &src = (*data(this->value_index(i)))[level];
+        const MGLevelObject<VECTOR>* src = data.read_ptr<MGLevelObject<VECTOR> >(this->value_index(i));
         VectorSlice<std::vector<std::vector<double> > > dst(values[i], component, n_comp);
-        fe.get_function_values(src, make_slice(index, start, size), dst, true);
+        fe.get_function_values((*src)[level], make_slice(index, start, size), dst, true);
       }
 
     for (unsigned int i=0; i<this->n_gradients(); ++i)
       {
-        const VECTOR &src = (*data(this->value_index(i)))[level];
+        const MGLevelObject<VECTOR>* src = data.read_ptr<MGLevelObject<VECTOR> >(this->value_index(i));
         VectorSlice<std::vector<std::vector<Tensor<1,dim> > > > dst(gradients[i], component, n_comp);
-        fe.get_function_gradients(src, make_slice(index, start, size), dst, true);
+        fe.get_function_gradients((*src)[level], make_slice(index, start, size), dst, true);
       }
 
     for (unsigned int i=0; i<this->n_hessians(); ++i)
       {
-        const VECTOR &src = (*data(this->value_index(i)))[level];
+        const MGLevelObject<VECTOR>* src = data.read_ptr<MGLevelObject<VECTOR> >(this->value_index(i));
         VectorSlice<std::vector<std::vector<Tensor<2,dim> > > > dst(hessians[i], component, n_comp);
-        fe.get_function_hessians(src, make_slice(index, start, size), dst, true);
+        fe.get_function_hessians((*src)[level], make_slice(index, start, size), dst, true);
       }
   }
 
@@ -237,7 +248,7 @@ namespace MeshWorker
   MGVectorData<VECTOR, dim, spacedim>::memory_consumption () const
   {
     std::size_t mem = VectorSelector::memory_consumption();
-    mem += sizeof (data);
+    mem += sizeof (this->data);
     return mem;
   }
 }
diff --git a/tests/algorithms/theta_01.cc b/tests/algorithms/theta_01.cc
new file mode 100644
index 0000000000..bb18059d1b
--- /dev/null
+++ b/tests/algorithms/theta_01.cc
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+// $Id: theta_timestepping.cc 30050 2013-07-18 20:03:31Z maier $
+//
+// Copyright (C) 2005 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// See documentation of ThetaTimestepping for documentation of this example
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+
+#include <deal.II/algorithms/operator.h>
+#include <deal.II/algorithms/theta_timestepping.h>
+
+#include <iostream>
+
+using namespace dealii;
+using namespace Algorithms;
+
+
+class Explicit
+  : public Operator<Vector<double> >
+{
+public:
+  Explicit(const FullMatrix<double> &matrix);
+  void operator() (AnyData &out, const AnyData &in);
+private:
+  SmartPointer<const FullMatrix<double>, Explicit> matrix;
+  FullMatrix<double> m;
+};
+
+
+class Implicit
+  : public Operator<Vector<double> >
+{
+public:
+  Implicit(const FullMatrix<double> &matrix);
+  void operator() (AnyData &out, const AnyData &in);
+private:
+  SmartPointer<const FullMatrix<double>, Implicit> matrix;
+  FullMatrix<double> m;
+};
+
+// End of declarations
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  FullMatrix<double> matrix(2);
+  matrix(0,0) = 0.;
+  matrix(1,1) = 0.;
+  matrix(0,1) = numbers::PI;
+  matrix(1,0) = -numbers::PI;
+
+  OutputOperator<Vector<double> > out;
+  out.initialize_stream(logfile);
+
+  Explicit op_explicit(matrix);
+  Implicit op_implicit(matrix);
+  ThetaTimestepping<Vector<double> > solver(op_explicit, op_implicit);
+  solver.timestep_control().start_step(0.1);
+  //solver.set_output(out);
+
+  Vector<double> value(2);
+  value(0) = 1.;
+  AnyData indata;
+  AnyData outdata;
+  Vector<double> *p = &value;
+  outdata.add(p, "value");
+  deallog << "Initial: " << value(0) << ' ' << value(1) << std::endl;
+  solver.notify(Events::initial);
+  solver(outdata, indata);
+  deallog << "Result: " << value(0) << ' ' << value(1)
+	  << " Norm " << value.l2_norm() << std::endl;
+}
+
+
+Explicit::Explicit(const FullMatrix<double> &M)
+  :
+  matrix(&M)
+{
+  m.reinit(M.m(), M.n());
+}
+
+
+void
+Explicit::operator() (AnyData &out, const AnyData &in)
+{
+  const double* step = in.read_ptr<double>("Timestep");
+  
+  if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size))
+    {
+      m.equ(-*step, *matrix);
+      for (unsigned int i=0; i<m.m(); ++i)
+        m(i,i) += 1.;
+    }
+  this->notifications.clear();
+  unsigned int i = in.find("Previous iterate");
+  m.vmult(*out.entry<Vector<double>*>(0), *in.read_ptr<Vector<double> >(i));
+}
+
+
+Implicit::Implicit(const FullMatrix<double> &M)
+  :
+  matrix(&M)
+{
+  m.reinit(M.m(), M.n());
+}
+
+
+void
+Implicit::operator() (AnyData &out, const AnyData &in)
+{
+  const double * step = in.read_ptr<double>("Timestep");
+  
+  if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size))
+    {
+      m.equ(*step, *matrix);
+      for (unsigned int i=0; i<m.m(); ++i)
+        m(i,i) += 1.;
+      m.gauss_jordan();
+    }
+  this->notifications.clear();
+
+  unsigned int i = in.find("Previous time");
+  m.vmult(*out.entry<Vector<double>*>(0), *in.read_ptr<Vector<double> >(i));
+}
+
+
diff --git a/tests/algorithms/theta_01.output b/tests/algorithms/theta_01.output
new file mode 100644
index 0000000000..7c4d9ddfc7
--- /dev/null
+++ b/tests/algorithms/theta_01.output
@@ -0,0 +1,13 @@
+
+DEAL::Initial: 1.00000 0
+DEAL:Theta::Time step:0.100000
+DEAL:Theta::Time step:0.200000
+DEAL:Theta::Time step:0.300000
+DEAL:Theta::Time step:0.400000
+DEAL:Theta::Time step:0.500000
+DEAL:Theta::Time step:0.600000
+DEAL:Theta::Time step:0.700000
+DEAL:Theta::Time step:0.800000
+DEAL:Theta::Time step:0.900000
+DEAL:Theta::Time step:1.00000
+DEAL::Result: -0.999676 0.0254599 Norm 1.00000