]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Modified how the PointValueHistory class supports vector fields to allow a more consi...
authorrapson <rapson@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 12 Jul 2012 21:57:18 +0000 (21:57 +0000)
committerrapson <rapson@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 12 Jul 2012 21:57:18 +0000 (21:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@25695 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/numerics/point_value_history.h
deal.II/source/numerics/point_value_history.cc
deal.II/source/numerics/point_value_history.inst.in

index 5d3b023cfa8571fd434153f6dd349d24c72dc786..fc90377f7c761f2f312d13b5855b586e739635e3 100644 (file)
@@ -101,9 +101,10 @@ namespace internal
  * just once at the start. However, this method is not suitable for
  * FiniteElement objects that do not assign dofs to actual mesh locations
  * (i.e. FEs without @ref GlossSupport "support points") or if adaptive mesh refinement is used.
- * The class will throw an exception if the finite element does not have support points
- * or any change in the @p dof_handler is identified. (The test is simply
- * whether the number of dofs changes.)
+ * The class will throw an exception if any change to the triangulation is made. Although the 
+ * listeners that detect this change could be used to recompute the nearest support points, 
+ * the location of the support point will most likely change slightly, making the interpretation
+ * of the data difficult, hence this is not implemented currently.
  *
  * <li> Secondly, @p evaluate_field_at_requested_location calls @p
  * VectorTools::point_value to compute values at the specific point requested.
@@ -119,16 +120,21 @@ namespace internal
  * code using adaptive refinement, but as the output will be meaningful (in
  * the sense that the quadrature point selected is guaranteed to remain in the
  * same vicinity, the class does not prevent the use of this method in adaptive
- * codes.
+ * codes. The class provides warnings in the output files if the mesh has changed.
  * </ol>
  *
- * When recording a new mnemonic name, the user must chose whether it should be
- * a strictly scalar field (@p add_scalar_field_name) or whether it should
- * depend on the @p FE (@p add_field_name). Calling @p add_field_name with a
- * @p dof_handler using a vector @p FE will result in a vector field. The vector
- * field will be of dimention @p FE.n_components(). Therefore, in this class the
- * distinction between vector and scalar is as related to the @p FE and not
- * whether a quantity is a logical vector.
+ * When recording a new mnemonic name, the user must supply a std::vector<unsigned int>
+ * components that lists the (generally FE) components to be extracted from the 
+ * given input. If the user simply wants to extract all the components, or the first
+ * n components then default parameters and overloaded functions mean that 
+ * the std::vector need not be explicitly supplied to the @p add_field_name method. 
+ * If the @p evaluate_field with a @p DataPostprocessor object is used, the components
+ * are interpreted as the component of the @p DataPostprocessor return vector. This 
+ * vector can be larger than the FE space.
+ *
+ * The class automatically generates names for the data stored based on the mnemoics 
+ * supplied. The methods @p add_component_names and @p add_independent_names allow 
+ * the user to provide lists of names to use instead if desired.
  *
  * Following is a little code snippet that shows a common usage of this class:
  *
@@ -184,24 +190,20 @@ class PointValueHistory
 
                                      /**
                                       * Constructor linking the class to a
-                                      * specific @p DoFHandler. Since this
+                                      * specific @p DoFHandler. This
                                       * class reads specific data from the @p
                                       * DoFHandler and stores it internally
                                       * for quick access (in particular dof
                                       * indices of closest neighbors to
                                       * requested points) the class is fairly
                                       * intolerant to changes to the @p
-                                      * DoFHandler.  Mesh refinement and @p
+                                      * DoFHandler if data at support points
+                                      * is required. Mesh refinement and @p
                                       * DoFRenumbering methods should be
                                       * performed before the @p add_points
                                       * method is called and adaptive grid
-                                      * refinement is not supported (unless it
-                                      * is completed before points are added).
-                                      * Changes to the @p DoFHandler are
-                                      * tested for by looking for a change in
-                                      * the number of dofs, which may not
-                                      * always occur so the user must be aware
-                                      * of these limitations.
+                                      * refinement is only supported by some
+                                      * methods.
                                       *
                                       * The user can store extra variables
                                       * which do not relate to mesh location
@@ -217,9 +219,7 @@ class PointValueHistory
                        const unsigned int n_independent_variables = 0);
 
                                      /**
-                                      * Copy constructor explicitly provided,
-                                      * although default copy operator would
-                                      * be sufficient.  This constructor can
+                                      * Copy constructor. This constructor can
                                       * be safely called with a @p
                                       * PointValueHistory object that contains
                                       * data, but this could be expensive and
@@ -228,9 +228,7 @@ class PointValueHistory
     PointValueHistory (const PointValueHistory & point_value_history);
 
                                      /**
-                                      * Assignment operator explicitly
-                                      * provided, although default assignment
-                                      * operator would be sufficient.  This
+                                      * Assignment operator. This
                                       * assignment operator can be safely
                                       * called once the class is closed and
                                       * data added, but this is provided
@@ -252,7 +250,7 @@ class PointValueHistory
                                       * Add a single point to the class. The
                                       * support points (one per component) in
                                       * the mesh that are closest to that
-                                      * point is found and its details stored
+                                      * point are found and their details stored
                                       * for use when @p evaluate_field is
                                       * called. If more than one point is
                                       * required rather used the @p add_points
@@ -265,7 +263,7 @@ class PointValueHistory
                                       * Add multiple points to the class. The
                                       * support points (one per component) in
                                       * the mesh that are closest to that
-                                      * point is found and its details stored
+                                      * point is found and their details stored
                                       * for use when @p evaluate_field is
                                       * called. If more than one point is
                                       * required, rather call this method as it
@@ -286,41 +284,59 @@ class PointValueHistory
                                      /**
                                       * Put another mnemonic string (and hence
                                       * @p VECTOR) into the class. This method
-                                      * adds storage space for FE.n_components()
-                                      * variables and should be used with methods
-                                      * that extract all FE components from a vector.
+                                      * adds storage space for components.size()
+                                      * variables. 
                                       * This also adds extra entries for points
                                       * that are already in the class, so
                                       * @p add_field_name and @p add_points can
                                       * be called in any order.
                                       */
-    void add_field_name(const std::string &vector_name);
-
+    void add_field_name(const std::string &vector_name, 
+                        const std::vector <unsigned int> &components = std::vector <unsigned int>());
+                        
                                      /**
                                       * Put another mnemonic string (and hence
                                       * @p VECTOR) into the class. This method
-                                      * only adds storage space for a single
-                                      * variable whether the FE is scalar or not
-                                      * so it can be used with
-                                      * @p DataPostprocessor objects that compute
-                                      * a scalar value(s) from a vector field.
+                                      * adds storage space for components.size()
+                                      * variables. 
                                       * This also adds extra entries for points
                                       * that are already in the class, so
                                       * @p add_field_name and @p add_points can
                                       * be called in any order.
-                                      */
-    void add_scalar_field_name(const std::string &vector_name);
+                                      * This method generates a std::vector 
+                                      * 0, ..., n_components-1 and calls the 
+                                      * previous function.
+                                      */
+    void add_field_name(const std::string &vector_name, 
+                        const unsigned int n_components);
+                        
+                                       /**
+                                        * Provide optional names for each component
+                                        * of a field. These names will be used 
+                                        * instead of names generated from the 
+                                        * field name, if supplied.
+                                        */
+    void add_component_names(const std::string &vector_name, 
+                        const std::vector <std::string> &component_names);
+                        
+                                       /**
+                                        * Provide optional names for the
+                                        * independent values. These names will 
+                                        * be used instead of "Indep_...", if
+                                        * supplied.
+                                        */
+    void add_independent_names(const std::vector <std::string> &independent_names);
+
+
 
                                      /**
                                       * Extract values at the stored points
                                       * from the VECTOR supplied and add them
                                       * to the new dataset in vector_name.
-                                      * The scalar_component input is only
-                                      * used if the FE has multiple components
-                                      * and the field to be evaluated is scalar.
-                                      * In this case, the scalar_component input
-                                      * is used to select the component used. If
-                                      * a @p DoFHandler is used, one (and only
+                                      * The list of components supplied when the
+                                      * field was added is used to select
+                                      * components to extract.
+                                      * If a @p DoFHandler is used, one (and only
                                       * one) evaluate_field method
                                       * must be called for each dataset (time
                                       * step, iteration, etc) for each
@@ -328,23 +344,30 @@ class PointValueHistory
                                       * ExcDataLostSync error can occur.
                                       */
     template <class VECTOR>
-    void evaluate_field(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0);
+    void evaluate_field(const std::string &vector_name, const VECTOR & solution);
 
 
                                      /**
-                                      * Compute values using a DataPostprocessor object
-                                      * with the VECTOR supplied and add them
+                                      * Compute values using a @p DataPostprocessor object
+                                      * with the @p VECTOR supplied and add them
                                       * to the new dataset in vector_name.
-                                      * The DataPostprocessor component_interpretation
-                                      * is used to determine select scalar or vector data.
-                                      * If the component_is_part_of_vector flag is used, then
-                                      * the indicated vector must have the same number of
-                                      * components as the FE associated with the dof_handler.
-                                      * The DataPostprocessor object is queried for field names.
-                                      * In the case of vector fields, the name of the first
-                                      * component is used as the overall name of the vector.
+                                      * The list of components supplied when the
+                                      * field was added is used to select
+                                      * components to extract from the @p DataPostprocessor
+                                      * return vector.
+                                      * This method takes a vector of field names to
+                                      * process and is preferred if many fields use the
+                                      * same @p DataPostprocessor object as each cell
+                                      * is only located once.
                                       * The quadrature object supplied is used for all
                                       * components of a vector field.
+                                      * Although this method will not throw an exception
+                                      * if the mesh has changed. (No internal data structures
+                                      * are invalidated as the quadrature points are repicked
+                                      * each time the function is called.) Nevertheless
+                                      * the user must be aware that if the mesh changes the
+                                      * point selected will also vary slightly, making 
+                                      * interpretation of the data more difficult. 
                                       * If a @p DoFHandler is used, one (and only
                                       * one) evaluate_field method
                                       * must be called for each dataset (time
@@ -353,7 +376,17 @@ class PointValueHistory
                                       * ExcDataLostSync error can occur.
                                       */
     template <class VECTOR>
-    void evaluate_field(const VECTOR & solution, const DataPostprocessor<dim> & data_postprocessor, const Quadrature<dim> & quadrature);
+    void evaluate_field(const std::vector <std::string> &vector_names, const VECTOR & solution, const DataPostprocessor<dim> & data_postprocessor, const Quadrature<dim> & quadrature);
+
+                                     /** 
+                                      * Construct a std::vector <std::string> 
+                                      * containing only vector_name and
+                                      * call the above function. The above function
+                                      * is more efficient if multiple fields
+                                      * use the same @p DataPostprocessor object.
+                                      */
+    template <class VECTOR>
+    void evaluate_field(const std::string &vector_name, const VECTOR & solution, const DataPostprocessor<dim> & data_postprocessor, const Quadrature<dim> & quadrature);
 
 
                                      /**
@@ -365,15 +398,11 @@ class PointValueHistory
                                       * has been modified because it uses calls to
                                       * @p VectorTools::point_value to extract there
                                       * data. Therefore, if only this method is used,
-                                      * the class can be used with adaptive refinement.
-                                      * Specifically, the vector form of
-                                      * @p VectorTools::point_value that assumes a
-                                      * Q1 mapping is used.
-                                      * The scalar_component input is only
-                                      * used if the FE has multiple components
-                                      * and the field to be evaluated is scalar.
-                                      * In this case, the scalar_component input
-                                      * is used to select the component used. If
+                                      * the class is fully compatible with 
+                                      * adaptive refinement.
+                                      * The list of components supplied when the
+                                      * field was added is used to select
+                                      * components to extract. If
                                       * a @p DoFHandler is used, one (and only
                                       * one) evaluate_field method
                                       * must be called for each dataset (time
@@ -382,7 +411,8 @@ class PointValueHistory
                                       * ExcDataLostSync error can occur.
                                       */
     template <class VECTOR>
-    void evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0);
+    void evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution); // changed
+
 
                                      /**
                                       * Add the key for the current dataset to
@@ -423,6 +453,10 @@ class PointValueHistory
                                       * key and independent data. The file
                                       * name agrees with the order the points
                                       * were added to the class.
+                                      * The names of the data columns can
+                                      * be supplied using the functions
+                                      * @p add_component_names and 
+                                      * @p add_independent_names.
                                       * The support point information is only
                                       * meaningful if the dof_handler has not
                                       * been changed. Therefore, if adaptive
@@ -527,7 +561,7 @@ class PointValueHistory
                                       * to the correct number of points by the
                                       * method.
                                       */
-    void get_postprocessor_locations (std::vector<Point <dim> > & locations, const Quadrature<dim> & quadrature);
+    void get_postprocessor_locations (const Quadrature<dim> & quadrature, std::vector<Point <dim> > & locations);
 
                                      /**
                                       * Once datasets have been added to the
@@ -593,6 +627,19 @@ class PointValueHistory
 
     bool deep_check (bool strict);
 
+
+                                       /**
+                                        *  A function that will be triggered
+                                        *  through signals whenever the
+                                        *  triangulation is modified.
+                                        *
+                                        *  It is currently used to check
+                                        *  if the triangulation has changed,
+                                        *  invalidating precomputed values.
+                                        */
+    void tria_change_listener ();
+
+
                                      /**
                                       * A call has been made to @p
                                       * push_back_independent when no
@@ -620,26 +667,21 @@ class PointValueHistory
                                      /**
                                       * A method which requires access to a @p
                                       * DoFHandler to be meaningful has been
-                                      * called when @p n_dofs is reported to
-                                      * be zero (most likely due to default
-                                      * constructor being called).  Only
+                                      * called when @p have_dof_handler is 
+                                      * false (most likely due to default
+                                      * constructor being called). Only
                                       * independent variables may be logged
                                       * with no DoFHandler.
                                       */
     DeclException0(ExcDoFHandlerRequired);
 
                                      /**
-                                      * @p n_dofs for the @p DoFHandler has
-                                      * changed since the constructor was
-                                      * called. Possibly the mesh has been
-                                      * refined in some way. This suggests
+                                      *    The triangulation indicated that mesh 
+                                      * has been refined in some way. This suggests
                                       * that the internal dof indices stored
                                       * are no longer meaningful.
                                       */
-    DeclException2(ExcDoFHandlerChanged,
-                   int,
-                   int,
-                   << "Original n_dofs (" << arg1 << ") != current n_dofs (" << arg2 << ")");
+    DeclException0(ExcDoFHandlerChanged);
 
   private:
                                      /**
@@ -654,6 +696,14 @@ class PointValueHistory
                                       * location.
                                       */
     std::vector <std::vector <double> > independent_values;
+    
+                                     /**
+                                      * Saves a vector listing component
+                                      * names associated with a independent_values.
+                                      * This will be an empty vector
+                                      * if the user does not supplies names.
+                                      */
+    std::vector<std::string> indep_names;
 
                                      /**
                                       * Saves data for each mnemonic entry.
@@ -667,14 +717,22 @@ class PointValueHistory
                                       * FE.n_components () long.
                                       */
     std::map <std::string, std::vector <std::vector <double> > > data_store;
-
+    
+    
                                      /**
-                                      * Saves a boolean representing whether
-                                      * the mnemonic links to a scalar field
-                                      * (true) or not (false).
+                                      * Saves a vector listing components
+                                      * associated with a mnemonic. 
                                       */
-    std::map <std::string, bool> scalar_field;
-
+    std::map <std::string, std::vector<unsigned int> > field_components;
+    
+                                     /**
+                                      * Saves a vector listing component
+                                      * names associated with a mnemonic.
+                                      * This will be an empty vector
+                                      * if the user does not supplies names.
+                                      */
+    std::map <std::string, std::vector<std::string> > component_names_map;
+    
                                      /**
                                       * Saves the location and other mesh info
                                       * about support points.
@@ -682,11 +740,6 @@ class PointValueHistory
     std::vector <internal::PointValueHistory::PointGeometryData <dim> >
     point_geometry_data;
 
-                                     /**
-                                      * A tempory pair used when adding to
-                                      * @p data_store.
-                                      */
-    std::pair<std::string, std::vector <std::vector <double> > > pair_data; // could possibly be removed
 
                                      /**
                                       * Used to enforce @p closed state for some
@@ -708,17 +761,26 @@ class PointValueHistory
                                       */
     SmartPointer<const DoFHandler<dim>,PointValueHistory<dim> > dof_handler;
 
+
                                      /**
-                                      * Variable to check that number of dofs
-                                      * in the dof_handler does not change.
-                                      *
-                                     * A cheap way to check that the
-                                      * triangulation has not been refined in
-                                      * anyway since the class was set up.
-                                      * Refinement will invalidate stored dof
-                                      * indices linked to support points.
+                                      * Variable to check if the triangulation
+                                      * has changed. If it has changed, certain
+                                      * data is out of date (especially the
+                                      * PointGeometryData::solution_indices.
+                                      */
+    bool triangulation_changed;
+
+                                     /**
+                                      * A boolean to record whether the class was
+                                      * initialized with a DoFHandler or not.
                                       */
-    unsigned int n_dofs;
+    bool have_dof_handler;
+
+                                     /**
+                                      * Used to detect signals from the Triangulation.
+                                      */
+    boost::signals2::connection tria_listener;
+
                                      /**
                                       * Stores the number of independent
                                       * variables requested.
index 2cfa0ea914ed8cf6354da78b9cb40ecf2d0a2403..208cdaabe99512f783552a6fa3035b33075f491d 100644 (file)
@@ -36,10 +36,10 @@ namespace internal
     template <int dim>
     PointGeometryData<dim>
     ::PointGeometryData (const Point <dim> &new_requested_location, 
-                        const std::vector <Point <dim> > &new_locations,
-                        const std::vector <int> &new_sol_indices)
+             const std::vector <Point <dim> > &new_locations,
+             const std::vector <int> &new_sol_indices)
     {
-         requested_location = new_requested_location;
+      requested_location = new_requested_location;
       support_point_locations = new_locations;
       solution_indices = new_sol_indices;
     }
@@ -51,36 +51,45 @@ namespace internal
 template <int dim>
 PointValueHistory<dim>
 ::PointValueHistory (const unsigned int n_independent_variables) :
-               n_dofs (0),
-               n_indep (n_independent_variables)
+        n_indep (n_independent_variables)
 {
   closed = false;
   cleared = false;
-                                  // make a vector for keys
+  triangulation_changed = false;
+  have_dof_handler = false;
+
+                   // make a vector for keys
   dataset_key = std::vector <double> (); // initialize the std::vector
 
-                                  // make a vector of independent values
+                   // make a vector of independent values
   independent_values
     = std::vector<std::vector <double> > (n_indep, std::vector <double> (0));
+  indep_names = std::vector <std::string> ();
 }
 
 
 
 template <int dim>
 PointValueHistory<dim>::PointValueHistory (const DoFHandler<dim> & dof_handler,
-                    const unsigned int n_independent_variables) :
-               dof_handler (&dof_handler),
-               n_dofs (dof_handler.n_dofs ()),
-               n_indep (n_independent_variables)
+             const unsigned int n_independent_variables) :
+        dof_handler (&dof_handler),
+        n_indep (n_independent_variables)
 {
   closed = false;
   cleared = false;
-                                  // make a vector to store keys
+  triangulation_changed = false;
+  have_dof_handler = true;
+
+                   // make a vector to store keys
   dataset_key = std::vector <double> (); // initialize the std::vector
 
-                                  // make a vector for the independent values
+                   // make a vector for the independent values
   independent_values
     = std::vector<std::vector <double> > (n_indep, std::vector <double> (0));
+  indep_names = std::vector <std::string> ();
+
+  tria_listener = dof_handler.get_tria().signals.any_change.connect (std_cxx1x::bind (&PointValueHistory<dim>::tria_change_listener,
+                                std_cxx1x::ref(*this)));
 }
 
 
@@ -90,17 +99,28 @@ PointValueHistory<dim>::PointValueHistory (const PointValueHistory & point_value
 {
   dataset_key = point_value_history.dataset_key;
   independent_values = point_value_history.independent_values;
+  indep_names = point_value_history.indep_names;
   data_store = point_value_history.data_store;
-  scalar_field = point_value_history.scalar_field;
+  field_components = point_value_history.field_components;
+  component_names_map = point_value_history.component_names_map;
   point_geometry_data = point_value_history.point_geometry_data;
-  pair_data = point_value_history.pair_data;
+  
   closed = point_value_history.closed;
   cleared = point_value_history.cleared;
 
   dof_handler = point_value_history.dof_handler;
 
-  n_dofs = point_value_history.n_dofs;
+  triangulation_changed = point_value_history.triangulation_changed;
+  have_dof_handler = point_value_history.have_dof_handler;
   n_indep = point_value_history.n_indep;
+
+                           // What to do with tria_listener? 
+                           // Presume subscribe new instance?
+  if (have_dof_handler)
+  {
+    tria_listener = dof_handler->get_tria().signals.any_change.connect (std_cxx1x::bind     (&PointValueHistory<dim>::tria_change_listener,
+                                std_cxx1x::ref(*this)));
+  }
 }
 
 
@@ -111,17 +131,29 @@ PointValueHistory<dim>::operator= (const PointValueHistory & point_value_history
 {
   dataset_key = point_value_history.dataset_key;
   independent_values = point_value_history.independent_values;
+  indep_names = point_value_history.indep_names;
   data_store = point_value_history.data_store;
-  scalar_field = point_value_history.scalar_field;
+  field_components = point_value_history.field_components;
+  component_names_map = point_value_history.component_names_map;
   point_geometry_data = point_value_history.point_geometry_data;
-  pair_data = point_value_history.pair_data;
+
   closed = point_value_history.closed;
   cleared = point_value_history.cleared;
 
   dof_handler = point_value_history.dof_handler;
 
-  n_dofs = point_value_history.n_dofs;
+  triangulation_changed = point_value_history.triangulation_changed;
+  have_dof_handler = point_value_history.have_dof_handler;
   n_indep = point_value_history.n_indep;
+
+                           // What to do with tria_listener? 
+                           // Presume subscribe new instance?
+  if (have_dof_handler)
+  {
+    tria_listener = dof_handler->get_tria().signals.any_change.connect (std_cxx1x::bind     (&PointValueHistory<dim>::tria_change_listener,
+                                std_cxx1x::ref(*this)));
+  }
+
   return * this;
 }
 
@@ -131,6 +163,10 @@ template <int dim>
 PointValueHistory<dim>
 ::~PointValueHistory ()
 {
+  if (have_dof_handler)
+  {
+    tria_listener.disconnect ();
+  }
 }
 
 
@@ -139,49 +175,48 @@ template <int dim>
 void PointValueHistory<dim>
 ::add_point (const Point <dim> & location)
 {
-                                  // can't be closed to add additional points
-                                  // or vectors
+                   // can't be closed to add additional points
+                   // or vectors
   AssertThrow (!closed, ExcInvalidState ());
   AssertThrow (!cleared, ExcInvalidState ());
-  AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
-  AssertThrow (n_dofs == dof_handler->n_dofs (),
-              ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+  AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
 
-                                  // Implementation assumes that support
-                                  // points locations are dofs locations
+                   // Implementation assumes that support
+                   // points locations are dofs locations
   AssertThrow (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ());
 
-                                  // FEValues object to extract quadrature
-                                  // points from
+                   // FEValues object to extract quadrature
+                   // points from
   std::vector <Point <dim> >
     unit_support_points = dof_handler->get_fe ().get_unit_support_points ();
 
-                                  // While in general quadrature points seems
-                                  // to refer to Gauss quadrature points, in
-                                  // this case the quadrature points are
-                                  // forced to be the support points of the
-                                  // FE.
+                   // While in general quadrature points seems
+                   // to refer to Gauss quadrature points, in
+                   // this case the quadrature points are
+                   // forced to be the support points of the
+                   // FE.
   Quadrature<dim>
     support_point_quadrature (dof_handler->get_fe ().get_unit_support_points ());
   FEValues<dim> fe_values (dof_handler->get_fe (),
-                          support_point_quadrature,
-                          update_quadrature_points);
+               support_point_quadrature,
+               update_quadrature_points);
   unsigned int n_support_points
     = dof_handler->get_fe ().get_unit_support_points ().size ();
   unsigned int n_components
     = dof_handler->get_fe ().n_components ();
 
-                                  // set up a loop over all the cells in the
-                                  // DoFHandler
+                   // set up a loop over all the cells in the
+                   // DoFHandler
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler->begin_active ();
   typename DoFHandler<dim>::active_cell_iterator
     endc = dof_handler->end ();
 
-                                  // default values to be replaced as closer
-                                  // points are found however they need to be
-                                  // consistent in case they are actually
-                                  // chosen
+                   // default values to be replaced as closer
+                   // points are found however they need to be
+                   // consistent in case they are actually
+                   // chosen
   typename DoFHandler<dim>::active_cell_iterator current_cell = cell;
   std::vector <unsigned int> current_fe_index (n_components, 0); // need one index per component
   fe_values.reinit (cell);
@@ -189,32 +224,32 @@ void PointValueHistory<dim>
   for (unsigned int support_point = 0;
        support_point < n_support_points; support_point++)
     {
-                                      // setup valid data in the empty
-                                      // vectors
+                       // setup valid data in the empty
+                       // vectors
       unsigned int component
-       = dof_handler->get_fe ().system_to_component_index (support_point).first;
+    = dof_handler->get_fe ().system_to_component_index (support_point).first;
       current_points [component] = fe_values.quadrature_point (support_point);
       current_fe_index [component] = support_point;
     }
 
-                                  // check each cell to find a suitable
-                                  // support points
+                   // check each cell to find a suitable
+                   // support points
   for (; cell != endc; cell++)
     {
       fe_values.reinit (cell);
 
       for (unsigned int support_point = 0;
-          support_point < n_support_points; support_point++)
+       support_point < n_support_points; support_point++)
         {
           unsigned int component
-           = dof_handler->get_fe ().system_to_component_index (support_point).first;
+        = dof_handler->get_fe ().system_to_component_index (support_point).first;
           Point<dim> test_point
-           = fe_values.quadrature_point (support_point);
+        = fe_values.quadrature_point (support_point);
 
           if (location.distance (test_point) <
-             location.distance (current_points [component]))
+          location.distance (current_points [component]))
             {
-                                              // save the data
+                           // save the data
               current_points [component] = test_point;
               current_cell = cell;
               current_fe_index [component] = support_point;
@@ -227,30 +262,30 @@ void PointValueHistory<dim>
     local_dof_indices (dof_handler->get_fe ().dofs_per_cell);
   std::vector <int> new_solution_indices;
   current_cell->get_dof_indices (local_dof_indices);
-                                  // there is an implicit assumption here
-                                  // that all the closest support point to
-                                  // the requested point for all finite
-                                  // element components lie in the same cell.
-                                  // this could possibly be violated if
-                                  // components use different fe orders,
-                                  // requested points are on the edge or
-                                  // vertex of a cell and we are unlucky with
-                                  // floating point rounding. Worst case
-                                  // scenario however is that the point
-                                  // selected isn't the closest possible, it
-                                  // will still lie within one cell distance.
-                                  // calling
-                                  // GridTools::find_active_cell_around_point
-                                  // to obtain a cell to search may be an
-                                  // option for these methods, but currently
-                                  // the GridTools method does not cater for
-                                  // a vector of points, and does not seem to
-                                  // be intrinsicly faster than this method.
+                   // there is an implicit assumption here
+                   // that all the closest support point to
+                   // the requested point for all finite
+                   // element components lie in the same cell.
+                   // this could possibly be violated if
+                   // components use different fe orders,
+                   // requested points are on the edge or
+                   // vertex of a cell and we are unlucky with
+                   // floating point rounding. Worst case
+                   // scenario however is that the point
+                   // selected isn't the closest possible, it
+                   // will still lie within one cell distance.
+                   // calling
+                   // GridTools::find_active_cell_around_point
+                   // to obtain a cell to search may be an
+                   // option for these methods, but currently
+                   // the GridTools method does not cater for
+                   // a vector of points, and does not seem to
+                   // be intrinsicly faster than this method.
   for (unsigned int component = 0;
        component < dof_handler->get_fe ().n_components (); component++)
     {
       new_solution_indices
-       .push_back (local_dof_indices[current_fe_index [component]]);
+    .push_back (local_dof_indices[current_fe_index [component]]);
     }
 
   internal::PointValueHistory::PointGeometryData<dim>
@@ -261,11 +296,10 @@ void PointValueHistory<dim>
     data_store_begin = data_store.begin ();
   for (; data_store_begin != data_store.end (); data_store_begin++)
     {
-                                      // add an extra row to each vector
-                                      // entry
-         bool scalar = (scalar_field.find (data_store_begin->first))->second;
-      for (unsigned int component = 0;
-          component < (scalar ? 1 : dof_handler->get_fe ().n_components ()); component++)
+                       // add an extra row to each vector
+                       // entry
+      unsigned int n_components = (field_components.find (data_store_begin->first))->second.size();
+      for (unsigned int component = 0; component < n_components; component++)
         {
           data_store_begin->second.push_back (std::vector<double> (0));
         }
@@ -278,46 +312,48 @@ template <int dim>
 void PointValueHistory<dim>
 ::add_points (const std::vector <Point <dim> > & locations)
 {
-                                  // This algorithm adds points in the same
-                                  // order as they appear in the vector
-                                  // locations and users may depend on this
-                                  // so do not change order added!
+                   // This algorithm adds points in the same
+                   // order as they appear in the vector
+                   // locations and users may depend on this
+                   // so do not change order added!
 
-                                  // can't be closed to add additional points or vectors
+                   // can't be closed to add additional points or vectors
   AssertThrow (!closed, ExcInvalidState ());
   AssertThrow (!cleared, ExcInvalidState ());
-  AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
-  AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
-                                  // Implementation assumes that support
-                                  // points locations are dofs locations
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+  AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
+
+  
+                   // Implementation assumes that support
+                   // points locations are dofs locations
   AssertThrow (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ());
 
-                                  // FEValues object to extract quadrature
-                                  // points from
+                   // FEValues object to extract quadrature
+                   // points from
   std::vector <Point <dim> > unit_support_points = dof_handler->get_fe ().get_unit_support_points ();
 
-                                  // While in general quadrature points seems
-                                  // to refer to Gauss quadrature points, in
-                                  // this case the quadrature points are
-                                  // forced to be the support points of the
-                                  // FE.
+                   // While in general quadrature points seems
+                   // to refer to Gauss quadrature points, in
+                   // this case the quadrature points are
+                   // forced to be the support points of the
+                   // FE.
   Quadrature<dim> support_point_quadrature (dof_handler->get_fe ().get_unit_support_points ());
   FEValues<dim> fe_values (dof_handler->get_fe (), support_point_quadrature, update_quadrature_points);
   unsigned int n_support_points = dof_handler->get_fe ().get_unit_support_points ().size ();
   unsigned int n_components = dof_handler->get_fe ().n_components ();
 
-                                  // set up a loop over all the cells in the
-                                  // DoFHandler
+                   // set up a loop over all the cells in the
+                   // DoFHandler
   typename DoFHandler<dim>::active_cell_iterator cell = dof_handler->begin_active ();
   typename DoFHandler<dim>::active_cell_iterator endc = dof_handler->end ();
 
-                                  // default values to be replaced as closer
-                                  // points are found however they need to be
-                                  // consistent in case they are actually
-                                  // chosen vector <vector>s defined where
-                                  // previously single vectors were used
+                   // default values to be replaced as closer
+                   // points are found however they need to be
+                   // consistent in case they are actually
+                   // chosen vector <vector>s defined where
+                   // previously single vectors were used
 
-                                  // need to store one value per point per component
+                   // need to store one value per point per component
   std::vector <typename DoFHandler<dim>::active_cell_iterator > current_cell (locations.size (), cell);
 
   fe_values.reinit (cell);
@@ -325,8 +361,8 @@ void PointValueHistory<dim>
   std::vector <unsigned int> temp_fe_index (n_components, 0);
   for (unsigned int support_point = 0; support_point < n_support_points; support_point++)
     {
-                                      // setup valid data in the empty
-                                      // vectors
+                       // setup valid data in the empty
+                       // vectors
       unsigned int component = dof_handler->get_fe ().system_to_component_index (support_point).first;
       temp_points [component] = fe_values.quadrature_point (support_point);
       temp_fe_index [component] = support_point;
@@ -334,8 +370,8 @@ void PointValueHistory<dim>
   std::vector <std::vector <Point <dim> > > current_points (locations.size (), temp_points); // give a valid start point
   std::vector <std::vector <unsigned int> > current_fe_index (locations.size (), temp_fe_index);
 
-                                  // check each cell to find suitable support
-                                  // points
+                   // check each cell to find suitable support
+                   // points
   for (; cell != endc; cell++)
     {
       fe_values.reinit (cell);
@@ -348,7 +384,7 @@ void PointValueHistory<dim>
             {
               if (locations[point].distance (test_point) < locations[point].distance (current_points[point][component]))
                 {
-                                                  // save the data
+                   // save the data
                   current_points[point][component] = test_point;
                   current_cell[point] = cell;
                   current_fe_index[point][component] = support_point;
@@ -373,13 +409,13 @@ void PointValueHistory<dim>
       point_geometry_data.push_back (new_point_geometry_data);
 
       std::map <std::string, std::vector <std::vector <double> > >::iterator
-       data_store_begin = data_store.begin ();
+      data_store_begin = data_store.begin ();
       for (; data_store_begin != data_store.end (); data_store_begin++)
         {
-                                          // add an extra row to each vector
-                                          // entry
-                 bool scalar = (scalar_field.find (data_store_begin->first))->second;
-          for (unsigned int component = 0; component < (scalar ? 1 : dof_handler->get_fe ().n_components ()); component++)
+                       // add an extra row to each vector
+                       // entry
+          unsigned int n_components = (field_components.find (data_store_begin->first))->second.size();
+          for (unsigned int component = 0; component < n_components; component++)
             {
               data_store_begin->second.push_back (std::vector<double> (0));
             }
@@ -394,58 +430,84 @@ void PointValueHistory<dim>
 
 template <int dim>
 void PointValueHistory<dim>
-::add_field_name (const std::string &vector_name)
+::add_field_name (const std::string &vector_name, const std::vector <unsigned int> &components)
 {
-                                  // can't be closed to add additional points
-                                  // or vectors
+                   // can't be closed to add additional points
+                   // or vectors
   AssertThrow (!closed, ExcInvalidState ());
-  AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
   AssertThrow (!cleared, ExcInvalidState ());
-  AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+  AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
 
 
-                                  // make and add a new vector
-                                  // point_geometry_data.size() long
+                    // Make a field_components list, if not supplied
+  std::vector <unsigned int> temp_components = components;
+  if (temp_components.size() == 0)
+  {
+                    // Generate component list 
+                    // 0,..,fe.n_components -1
+      for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+          temp_components.push_back (component);
+  }
+  field_components.insert (std::pair <std::string, std::vector<unsigned int> > (vector_name, temp_components)); 
+  
+                   // insert an empty vector of strings
+                   // to ensure each field has an entry
+                   // in the map
+  std::pair <std::string, std::vector <std::string> > empty_names (vector_name, std::vector <std::string> ());
+  component_names_map.insert (empty_names); 
+  
+                   // make and add a new vector
+                   // point_geometry_data.size() long
+  std::pair<std::string, std::vector <std::vector <double> > > pair_data;
   pair_data.first = vector_name;
-  int n_datastreams = point_geometry_data.size () * (dof_handler->get_fe ().n_components ()); // each point has n_components sub parts
+  int n_datastreams = point_geometry_data.size () * temp_components.size(); // each point has field_components sub parts
   std::vector < std::vector <double> > vector_size (n_datastreams,
                                                     std::vector <double> (0));
   pair_data.second = vector_size;
   data_store.insert (pair_data);
-                                  // this is only a scalar field if 
-                                  // the FE is scalar
-  if (dof_handler->get_fe ().n_components () == 1)
-       scalar_field.insert (std::pair <std::string, bool> (vector_name, true));
-  else
-       scalar_field.insert (std::pair <std::string, bool> (vector_name, false));
 }
 
 
-
 template <int dim>
 void PointValueHistory<dim>
-::add_scalar_field_name (const std::string &vector_name)
+::add_field_name(const std::string &vector_name, const unsigned int n_components)
 {
-                                  // can't be closed to add additional points
-                                  // or vectors
-  AssertThrow (!closed, ExcInvalidState ());
-  AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
-  AssertThrow (!cleared, ExcInvalidState ());
-  AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+  std::vector <unsigned int> temp_components;
+                           // Generate component list 
+                           // 0,..,n_components -1
+  for (unsigned int component = 0; component < n_components; component++)
+    temp_components.push_back (component);
+  
+  add_field_name (vector_name, temp_components);
+}
 
 
-                                  // make and add a new vector
-                                  // point_geometry_data.size() long
-  pair_data.first = vector_name;
-  int n_datastreams = point_geometry_data.size (); // each point has 1 value for this function
-  std::vector < std::vector <double> > vector_size (n_datastreams,
-                                                    std::vector <double> (0));
-  pair_data.second = vector_size;
-  data_store.insert (pair_data);
-  scalar_field.insert (std::pair <std::string, bool> (vector_name, true)); // scalar by definition in this method
+template <int dim>
+void PointValueHistory<dim>
+::add_component_names(const std::string &vector_name, 
+                        const std::vector <std::string> &component_names)
+{
+  typename std::map <std::string, std::vector <std::string> >::iterator names = component_names_map.find(vector_name);
+  Assert (names != component_names_map.end(), ExcMessage("vector_name not in class"));
+
+  typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
+  Assert (components != field_components.end(), ExcMessage("vector_name not in class"));  
+  Assert (component_names.size() == components->second.size(), ExcDimensionMismatch (component_names.size(), components->second.size()));
+  
+  names->second = component_names;
 }
 
 
+template <int dim>
+void PointValueHistory<dim>
+::add_independent_names(const std::vector <std::string> &independent_names)
+{
+  Assert (independent_names.size() == n_indep, ExcDimensionMismatch (independent_names.size(), n_indep));
+  
+  indep_names = independent_names;
+}
+
 
 template <int dim>
 void PointValueHistory<dim>
@@ -462,6 +524,7 @@ void PointValueHistory<dim>
 {
   cleared = true;
   dof_handler = 0;
+  have_dof_handler = false;
 }
 
 // Need to test that the internal data has a full and complete dataset for
@@ -479,54 +542,49 @@ void PointValueHistory<dim>
 template <int dim>
 template <class VECTOR>
 void PointValueHistory<dim>
-::evaluate_field (const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component)
+::evaluate_field (const std::string &vector_name, const VECTOR & solution)
 {
-                                  // must be closed to add data to internal
-                                  // members.
+                           // must be closed to add data to internal
+                           // members.
   Assert (closed, ExcInvalidState ());
   Assert (!cleared, ExcInvalidState ());
-  Assert (n_dofs != 0, ExcDoFHandlerRequired ());
-  Assert (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+  AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
+  
   if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0
     {
       Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
     }
-                                  // Look up the field name and get an
-                                  // iterator for the map. Doing this
-                                  // up front means that it only needs
-                                  // to be done once and also allows us
-                                  // to check vector_name is in the map.
+                           // Look up the field name and get an
+                           // iterator for the map. Doing this
+                           // up front means that it only needs
+                           // to be done once and also allows us
+                           // to check vector_name is in the map.
   typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(vector_name);
   Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
+                   // Repeat for field_components
+  typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
+  Assert (components != field_components.end(), ExcMessage("vector_name not in class"));  
+  
+  Assert (components->second.size () <= dof_handler->get_fe ().n_components (), ExcIndexRange (components->second.size (), 0, dof_handler->get_fe ().n_components () + 1)); // 0 components legal as request to generate list internally!
+  for (unsigned int component = 0; component < components->second.size (); component++)
+  {
+      Assert (components->second [component] < dof_handler->get_fe ().n_components (), ExcIndexRange (components->second[component], 0, dof_handler->get_fe ().n_components ()));
+  } // hopefully the compiler will optimise 
+  
   typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
   for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
+  {
+                       // Look up the component to add 
+                       // in field_components, and
+                       // access the data associated with
+                       // that component
+    for (unsigned int comp = 0; comp < components->second.size(); comp++)
     {
-         if (!(scalar_field.find (vector_name))->second)
-         {
-                                      // The solution_indices data has a
-                                          // one-to-one correspondance to FE
-                                          // components therefore we simply
-                                      // step through each point, to access
-                                      // the Solution_index and the vector
-                                      // index
-        for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
-          {
-            unsigned int solution_index = point->solution_indices[component];
-            data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (solution (solution_index));
-          }
-         }
-         else 
-         {
-                                      // We have found a scalar field 
-                                          // added to a dof_handler with a
-                                          // vector FE. Therefore, we take
-                                          // the value specified by the user
-                                          // in the scalar_component input.
-                Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError());
-                unsigned int solution_index = point->solution_indices[scalar_component];
-         data_store_field->second[data_store_index].push_back (solution (solution_index));
-         }
+        unsigned int solution_index = point->solution_indices[components->second[comp]];
+        data_store_field->second[data_store_index * components->second.size() + comp].push_back (solution (solution_index));
     }
+  }
 }
 
 
@@ -536,243 +594,225 @@ void PointValueHistory<dim>
 template <int dim>
 template <class VECTOR>
 void PointValueHistory<dim>
-::evaluate_field(const VECTOR & solution, const DataPostprocessor< dim> & data_postprocessor, const Quadrature<dim> & quadrature)
+::evaluate_field(const std::vector <std::string> &vector_names, const VECTOR & solution, const DataPostprocessor< dim> & data_postprocessor, const Quadrature<dim> & quadrature)
 {
-                                  // must be closed to add data to internal
-                                  // members.
+                           // must be closed to add data to internal
+                           // members.
   Assert (closed, ExcInvalidState ());
   Assert (!cleared, ExcInvalidState ());
-  Assert (n_dofs != 0, ExcDoFHandlerRequired ());
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
   if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0
     {
       Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
     }
 
-// Make an FEValues object
+                           // Make an FEValues object
   const UpdateFlags update_flags = data_postprocessor.get_needed_update_flags() | update_quadrature_points;
   Assert (!(update_flags & update_normal_vectors),
-         ExcMessage("The update of normal vectors may not be requested for evaluation of "
-                    "data on cells via DataPostprocessor."));
+      ExcMessage("The update of normal vectors may not be requested for evaluation of "
+             "data on cells via DataPostprocessor."));
   FEValues<dim> fe_values (dof_handler->get_fe (), quadrature, update_flags);
   unsigned int n_components = dof_handler->get_fe ().n_components ();
   unsigned int n_quadrature_points = quadrature.size();
 
-   // extract the data names, used to size the data_postprocessor output
-   std::vector<std::string> field_names = data_postprocessor.get_names ();
+  unsigned int n_output_variables = data_postprocessor.get_names().size();
 
-// Loop over points and find correct cell
+  // Loop over points and find correct cell
   typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
   for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
     {
-               // we now have a point to query, need to know what cell it is in
-               Point <dim> requested_location = point->requested_location;
-               typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1<dim>(), *dof_handler, requested_location).first;
-
-
-               fe_values.reinit (cell);
-               std::vector< Vector< double > > computed_quantities (1, Vector <double> (field_names.size())); // just one point needed
-               
-               // The case of a scalar FE
-               if (n_components == 1)
-               {
-                 // Extract data for the PostProcessor object
-                 std::vector< double > uh (n_quadrature_points, 0.0);
-                 std::vector< Tensor< 1, dim > > duh (n_quadrature_points, Tensor <1, dim> ());
-                 std::vector< Tensor< 2, dim > > dduh (n_quadrature_points, Tensor <2, dim> ());
-                 std::vector<Point<dim> > dummy_normals (1, Point<dim> ());
-                 std::vector<Point<dim> > evaluation_points;
-                                                       // at each point there is
-                                                       // only one component of
-                                                       // value, gradient etc.
-                 if (update_flags & update_values)
-                       fe_values.get_function_values (solution,
-                                                                 uh);
-                 if (update_flags & update_gradients)
-                       fe_values.get_function_gradients (solution,
-                                                                         duh);
-                 if (update_flags & update_hessians)
-                       fe_values.get_function_hessians (solution,
-                                                                       dduh);
-
-                 // find the closest quadrature point
-                 evaluation_points = fe_values.get_quadrature_points();
-                 double distance = cell->diameter ();
-                 unsigned int selected_point = 0;
-                 for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
-                 {
-                         if (requested_location.distance (evaluation_points[q_point]) < distance)
-                         {
-                                 selected_point = q_point;
-                                 distance = requested_location.distance (evaluation_points[q_point]);
-                         }
-                 }
-
-                 // Call compute_derived_quantities_vector or compute_derived_quantities_scalar
-                 data_postprocessor.
-                       compute_derived_quantities_scalar(std::vector< double > (1, uh[selected_point]),
-                                                         std::vector< Tensor< 1, dim > > (1, duh[selected_point]),
-                                                         std::vector< Tensor< 2, dim > > (1, dduh[selected_point]),
-                                                         dummy_normals,
-                                                         std::vector<Point<dim> > (1, evaluation_points[selected_point]),
-                                                         computed_quantities);
-                                                          
-               }
-                 else  // The case of a vector FE
-               {
-                 // Extract data for the PostProcessor object
-                 std::vector< Vector< double > > uh (n_quadrature_points, Vector <double> (n_components));
-                 std::vector< std::vector< Tensor< 1, dim > > > duh (n_quadrature_points, std::vector< Tensor< 1, dim > > (n_components,  Tensor< 1, dim >()));
-                 std::vector< std::vector< Tensor< 2, dim > > > dduh (n_quadrature_points, std::vector< Tensor< 2, dim > > (n_components,  Tensor< 2, dim >()));
-                 std::vector<Point<dim> > dummy_normals  (1, Point<dim> ());
-                 std::vector<Point<dim> > evaluation_points;
-                                                       // at each point there is
-                                                       // a vector valued
-                                                       // function and its
-                                                       // derivative...
-                 if (update_flags & update_values)
-                       fe_values.get_function_values (solution,
-                                                                 uh);
-                 if (update_flags & update_gradients)
-                       fe_values.get_function_gradients (solution,
-                                                                       duh);
-                 if (update_flags & update_hessians)
-                       fe_values.get_function_hessians (solution,
-                                                                       dduh);
-
-                 // find the closest quadrature point
-                 evaluation_points = fe_values.get_quadrature_points();
-                 double distance = cell->diameter ();
-                 unsigned int selected_point = 0;
-                 for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
-                 {
-                         if (requested_location.distance (evaluation_points[q_point]) < distance)
-                         {
-                                 selected_point = q_point;
-                                 distance = requested_location.distance (evaluation_points[q_point]);
-                         }
-                 }
-
-                 // Call compute_derived_quantities_vector or compute_derived_quantities_scalar
-                 data_postprocessor.
-                       compute_derived_quantities_vector(std::vector< Vector< double > > (1, uh[selected_point]),
-                                                         std::vector< std::vector< Tensor< 1, dim > > > (1, duh[selected_point]),
-                                                         std::vector< std::vector< Tensor< 2, dim > > > (1, dduh[selected_point]),
-                                                         dummy_normals,
-                                                         std::vector<Point<dim> > (1, evaluation_points[selected_point]),
-                                                         computed_quantities);
-               }
-               
-         // we now have the data and need to save it
-//                     std::vector<DataComponentInterpretation::DataComponentInterpretation> component_interpretation = data_postprocessor.get_data_component_interpretation ();
-               
-               // loop over data names
-               typename std::vector<std::string>::iterator name = field_names.begin();
-               for (unsigned int computed_quantities_index = 0; name != field_names.end();)
-               {
-                                                               // this loop is incremented in
-                                                               // separate locations depending on
-                                                               // whether the name is a vector_field
-                                                               // scalar_field or not in the 
-                                                               // data_store at all.
-                       typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(*name);
-                       if (data_store_field != data_store.end()) // Ignore names not in the class
-                       {  
-                               if (!(scalar_field.find (*name))->second)
-                               {
-                                                               // The solution_indices data has a
-                                                               // one-to-one correspondance to FE
-                                                               // components therefore we simply
-                                                               // step through each point, to access
-                                                               // the Solution_index and the vector
-                                                               // index
-                                 for (unsigned int component = 0; component < n_components; component++, name++, computed_quantities_index++)
-                                       {
-                                         Assert (computed_quantities_index < field_names.size(), ExcMessage ("Insufficient components in DataPostprocessor for current vector field."));
-                                         
-                                         data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (computed_quantities[0](computed_quantities_index));
-                                       }
-                               }
-                               else 
-                               {
-                                                               // We have found a scalar field 
-                                                               // added to a dof_handler with a
-                                                               // vector FE. Therefore, we take
-                                                               // the value specified by the user
-                                                               // in the scalar_component input.
-                                 Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError());
-                                 data_store_field->second[data_store_index].push_back (computed_quantities[0](computed_quantities_index));
-                                 name ++;
-                                 computed_quantities_index ++;
-                               }
-                       }
-                       else
-                       {
-                               name ++;
-                               computed_quantities_index ++;
-                       }
-               }
+                           // we now have a point to query,
+                           // need to know what cell it is in
+        Point <dim> requested_location = point->requested_location;
+        typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1<dim>(), *dof_handler, requested_location).first;
+
+
+        fe_values.reinit (cell);
+        std::vector< Vector< double > > computed_quantities (1, Vector <double> (n_output_variables)); // just one point needed
+        
+                           // The case of a scalar FE
+        if (n_components == 1)
+        {
+                           // Extract data for the 
+                           // PostProcessor object
+          std::vector< double > uh (n_quadrature_points, 0.0);
+          std::vector< Tensor< 1, dim > > duh (n_quadrature_points, Tensor <1, dim> ());
+          std::vector< Tensor< 2, dim > > dduh (n_quadrature_points, Tensor <2, dim> ());
+          std::vector<Point<dim> > dummy_normals (1, Point<dim> ());
+          std::vector<Point<dim> > evaluation_points;
+                            // at each point there is
+                            // only one component of
+                            // value, gradient etc.
+          if (update_flags & update_values)
+            fe_values.get_function_values (solution,
+                                  uh);
+          if (update_flags & update_gradients)
+            fe_values.get_function_gradients (solution,
+                                      duh);
+          if (update_flags & update_hessians)
+            fe_values.get_function_hessians (solution,
+                                    dduh);
+
+                           // find the closest quadrature point
+          evaluation_points = fe_values.get_quadrature_points();
+          double distance = cell->diameter ();
+          unsigned int selected_point = 0;
+          for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
+          {
+              if (requested_location.distance (evaluation_points[q_point]) < distance)
+              {
+                  selected_point = q_point;
+                  distance = requested_location.distance (evaluation_points[q_point]);
+              }
+          }
+
+                           // Call compute_derived_quantities_vector 
+                           // or compute_derived_quantities_scalar
+          data_postprocessor.
+            compute_derived_quantities_scalar(std::vector< double > (1, uh[selected_point]),
+                              std::vector< Tensor< 1, dim > > (1, duh[selected_point]),
+                              std::vector< Tensor< 2, dim > > (1, dduh[selected_point]),
+                              dummy_normals,
+                              std::vector<Point<dim> > (1, evaluation_points[selected_point]),
+                              computed_quantities);
+                               
+        }
+          else     // The case of a vector FE
+        {
+                           // Extract data for the PostProcessor object
+          std::vector< Vector< double > > uh (n_quadrature_points, Vector <double> (n_components));
+          std::vector< std::vector< Tensor< 1, dim > > > duh (n_quadrature_points, std::vector< Tensor< 1, dim > > (n_components,  Tensor< 1, dim >()));
+          std::vector< std::vector< Tensor< 2, dim > > > dduh (n_quadrature_points, std::vector< Tensor< 2, dim > > (n_components,  Tensor< 2, dim >()));
+          std::vector<Point<dim> > dummy_normals  (1, Point<dim> ());
+          std::vector<Point<dim> > evaluation_points;
+                            // at each point there is
+                            // a vector valued
+                            // function and its
+                            // derivative...
+          if (update_flags & update_values)
+            fe_values.get_function_values (solution,
+                                  uh);
+          if (update_flags & update_gradients)
+            fe_values.get_function_gradients (solution,
+                                    duh);
+          if (update_flags & update_hessians)
+            fe_values.get_function_hessians (solution,
+                                    dduh);
+
+                           // find the closest quadrature point
+          evaluation_points = fe_values.get_quadrature_points();
+          double distance = cell->diameter ();
+          unsigned int selected_point = 0;
+          for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
+          {
+              if (requested_location.distance (evaluation_points[q_point]) < distance)
+              {
+                  selected_point = q_point;
+                  distance = requested_location.distance (evaluation_points[q_point]);
+              }
+          }
+
+                           // Call compute_derived_quantities_vector
+                           // or compute_derived_quantities_scalar
+          data_postprocessor.
+            compute_derived_quantities_vector(std::vector< Vector< double > > (1, uh[selected_point]),
+                              std::vector< std::vector< Tensor< 1, dim > > > (1, duh[selected_point]),
+                              std::vector< std::vector< Tensor< 2, dim > > > (1, dduh[selected_point]),
+                              dummy_normals,
+                              std::vector<Point<dim> > (1, evaluation_points[selected_point]),
+                              computed_quantities);
+        }
+        
+        
+                           // we now have the data and need to save it
+                           // loop over data names
+        typename std::vector<std::string>::const_iterator name = vector_names.begin();
+        for (; name != vector_names.end(); name++)
+        {
+            typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(*name);
+            Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
+                            // Repeat for field_components
+            typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(*name);
+            Assert (components != field_components.end(), ExcMessage("vector_name not in class")); 
+                                
+                            // Push back components.size() computed 
+                            // quantities. Each component requested
+                            // must be less than the postprocessor
+                            // n_output_variables.
+            for (unsigned int comp = 0; comp < components->second.size(); comp++)
+            {
+                Assert (components->second[comp] < n_output_variables, ExcIndexRange (components->second[comp], 0, n_output_variables));
+                data_store_field->second[data_store_index * components->second.size() + comp].push_back (computed_quantities[0](components->second[comp]));
+            }
+        }
     } // end of loop over points
 }
 
 
-
+template <int dim>
+template <class VECTOR>
+void PointValueHistory<dim>
+::evaluate_field(const std::string &vector_name, const VECTOR & solution, const DataPostprocessor<dim> & data_postprocessor, const Quadrature<dim> & quadrature)
+{
+    std::vector <std::string> vector_names;
+    vector_names.push_back (vector_name);
+    evaluate_field (vector_names, solution, data_postprocessor, quadrature);
+}
 
 
 
 template <int dim>
 template <class VECTOR>
 void PointValueHistory<dim>
-::evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component)
+::evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution)
 {
-                                  // must be closed to add data to internal
-                                  // members.
+                     // must be closed to add data to internal
+                   // members.
   Assert (closed, ExcInvalidState ());
   Assert (!cleared, ExcInvalidState ());
-  Assert (n_dofs != 0, ExcDoFHandlerRequired ());
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+  
   if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0
     {
       Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
     }
-                                  // Look up the field name and get an
-                                  // iterator for the map. Doing this
-                                  // up front means that it only needs
-                                  // to be done once and also allows us
-                                  // to check vector_name is in the map.
+                   // Look up the field name and get an
+                   // iterator for the map. Doing this
+                   // up front means that it only needs
+                   // to be done once and also allows us
+                   // to check vector_name is in the map.
   typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(vector_name);
   Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
+                   // Repeat for field_components
+  typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
+  Assert (components != field_components.end(), ExcMessage("vector_name not in class"));  
+  
+  Assert (components->second.size () <= dof_handler->get_fe ().n_components (), ExcIndexRange (components->second.size (), 0, dof_handler->get_fe ().n_components () + 1)); // 0 components legal as request to generate list internally!
+  for (unsigned int component = 0; component < components->second.size (); component++)
+  {
+      Assert (components->second [component] < dof_handler->get_fe ().n_components (), ExcIndexRange (components->second[component], 0, dof_handler->get_fe ().n_components ()));
+  } // hopefully the compiler will optimise 
+  
+  
   typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
   Vector <double> value (dof_handler->get_fe().n_components());
   for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
+  {
+                      // Make a Vector <double> for the value
+                      // at the point. It will have as many
+                      // components as there are in the fe.
+      VectorTools::point_value (*dof_handler, solution, point->requested_location, value);
+      
+          
+                       // Look up the component to add 
+                       // in field_components, and
+                       // access the data associated with
+                       // that component
+    for (unsigned int comp = 0; comp < components->second.size(); comp++)
     {
-                                         // Make a Vector <double> for the value
-                                         // at the point. It will have as many
-                                         // components as there are in the fe.
-         VectorTools::point_value (*dof_handler, solution, point->requested_location, value);
-         
-         if (!(scalar_field.find (vector_name))->second)
-         {
-                                      // The solution_indices data has a
-                                          // one-to-one correspondance to FE
-                                          // components therefore we simply
-                                      // step through each point, to access
-                                      // the Solution_index and the vector
-                                      // index
-        for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
-          {
-            data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (value (component));
-          }
-         }
-         else 
-         {
-                                      // We have found a scalar field 
-                                          // added to a dof_handler with a
-                                          // vector FE. Therefore, we take
-                                          // the value specified by the user
-                                          // in the scalar_component input.
-                Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError());
-         data_store_field->second[data_store_index].push_back (value (scalar_component));
-         }
+        data_store_field->second[data_store_index * components->second.size() + comp].push_back (value (components->second[comp]));
     }
+  }
 }
 
 
@@ -780,8 +820,8 @@ template <int dim>
 void PointValueHistory<dim>
 ::start_new_dataset (double key)
 {
-                                  // must be closed to add data to internal
-                                  // members.
+                   // must be closed to add data to internal
+                   // members.
   Assert (closed, ExcInvalidState ());
   Assert (!cleared, ExcInvalidState ());
   Assert (deep_check (false), ExcDataLostSync ());
@@ -795,8 +835,8 @@ template <int dim>
 void PointValueHistory<dim>
 ::push_back_independent (const std::vector <double> &indep_values)
 {
-                                  // must be closed to add data to internal
-                                  // members.
+                   // must be closed to add data to internal
+                   // members.
   Assert (closed, ExcInvalidState ());
   Assert (!cleared, ExcInvalidState ());
   Assert (indep_values.size () == n_indep, ExcDimensionMismatch (indep_values.size (), n_indep));
@@ -817,7 +857,7 @@ void PointValueHistory<dim>
   AssertThrow (!cleared, ExcInvalidState ());
   AssertThrow (deep_check (true), ExcDataLostSync ());
 
-                                  // write inputs to a file
+                   // write inputs to a file
   if (n_indep != 0)
     {
       std::string filename = base_name + "_indep.gpl";
@@ -825,16 +865,26 @@ void PointValueHistory<dim>
 
       to_gnuplot << "# Data independent of mesh location\n";
 
-                                      // write column headings
+                       // write column headings
       to_gnuplot << "# <Key> ";
 
-      for (unsigned int component = 0; component < n_indep; component++)
-        {
-          to_gnuplot << "<Indep_" << component << "> ";
-        }
-      to_gnuplot << "\n";
-
-                                      // write general data stored
+      if (indep_names.size() > 0)
+      {
+          for (unsigned int name = 0; name < indep_names.size(); name++)
+          {
+              to_gnuplot << "<" << indep_names [name] << "> ";
+          }
+          to_gnuplot << "\n";
+      }
+      else
+      {    
+          for (unsigned int component = 0; component < n_indep; component++)
+          {
+              to_gnuplot << "<Indep_" << component << "> ";
+          }
+          to_gnuplot << "\n";
+      }
+                       // write general data stored
       for (unsigned int key = 0; key < dataset_key.size (); key++)
         {
           to_gnuplot << dataset_key[key];
@@ -851,81 +901,104 @@ void PointValueHistory<dim>
 
 
 
-                                  // write points to a file
-  if (n_dofs != 0)
+                   // write points to a file
+  if (have_dof_handler)
     {
-      AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
-         AssertThrow (postprocessor_locations.size() == 0 || postprocessor_locations.size() == point_geometry_data.size(), ExcDimensionMismatch (postprocessor_locations.size(), point_geometry_data.size()));
-                                         // We previously required the 
-                                         // number of dofs to remain the
-                                         // same to provide some sort of
-                                         // test on the relevance of the 
-                                         // support point indices stored.
-                                         // We now relax that to allow 
-                                         // adaptive refinement strategies
-                                         // to make use of the
-                                         // evaluate_field_requested_locations
-                                         // method. Note that the support point
-                                         // information is not meaningful if
-                                         // the number of dofs has changed.
-      //AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
-
-      unsigned int n_components = dof_handler->get_fe ().n_components ();
+        AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+      AssertThrow (postprocessor_locations.size() == 0 || postprocessor_locations.size() == point_geometry_data.size(), ExcDimensionMismatch (postprocessor_locations.size(), point_geometry_data.size()));
+                      // We previously required the 
+                      // number of dofs to remain the
+                      // same to provide some sort of
+                      // test on the relevance of the 
+                      // support point indices stored.
+                      // We now relax that to allow 
+                      // adaptive refinement strategies
+                      // to make use of the
+                      // evaluate_field_requested_locations
+                      // method. Note that the support point
+                      // information is not meaningful if
+                      // the number of dofs has changed.
+      //AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
 
       typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
       for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
         {
-                                          // for each point, open a file to
-                                          // be written to
+                       // for each point, open a file to
+                       // be written to
           std::string filename = base_name + "_" + Utilities::int_to_string (data_store_index, 2) + ".gpl"; // store by order pushed back
-                                          // due to
-                                          // Utilities::int_to_string(data_store_index,
-                                          // 2) call, can handle up to 100
-                                          // points
+                       // due to
+                       // Utilities::int_to_string(data_store_index,
+                       // 2) call, can handle up to 100
+                       // points
           std::ofstream to_gnuplot (filename.c_str ());
 
-                                          // put helpful info about the
-                                          // support point into the file as
-                                          // comments
-                 to_gnuplot << "# Requested location: " << point->requested_location << "\n";
+                       // put helpful info about the
+                       // support point into the file as
+                       // comments
+          to_gnuplot << "# Requested location: " << point->requested_location << "\n";
           to_gnuplot << "# DoF_index : Support location (for each component)\n";
-          for (unsigned int component = 0; component < n_components; component++)
+          for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
             {
               to_gnuplot << "# " << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n";
             }
-                 if (postprocessor_locations.size() != 0)
-                 {
-                         to_gnuplot << "# Postprocessor location: " << postprocessor_locations[data_store_index] << "\n";
-                 }
+          if (triangulation_changed)
+              to_gnuplot << "# (Original components and locations, may be invalidated by mesh change.)\n";
+
+          if (postprocessor_locations.size() != 0)
+          {
+              to_gnuplot << "# Postprocessor location: " << postprocessor_locations[data_store_index];
+              if (triangulation_changed)
+                  to_gnuplot << " (may be approximate)\n";
+          }
           to_gnuplot << "#\n";
 
-                                          // write column headings
+          
+                       // write column headings
           std::map <std::string, std::vector <std::vector <double> > >::iterator
-                         data_store_begin = data_store.begin ();
+              data_store_begin = data_store.begin ();
           to_gnuplot << "# <Key> ";
 
-          for (unsigned int component = 0; component < n_indep; component++)
-            {
-              to_gnuplot << "<Indep_" << component << "> ";
-            }
+          if (indep_names.size() > 0)
+          {
+              for (unsigned int name = 0; name < indep_names.size(); name++)
+              {
+                  to_gnuplot << "<" << indep_names [name] << "> ";
+              }
+          }
+          else
+          {    
+              for (unsigned int component = 0; component < n_indep; component++)
+              {
+                  to_gnuplot << "<Indep_" << component << "> ";
+              }
+          }
 
           for (; data_store_begin != data_store.end (); data_store_begin++)
             {
-                         if (!(scalar_field.find (data_store_begin->first))->second)
-                for (unsigned int component = 0; component < n_components; component++)
-                  {
-                    to_gnuplot << "<" << data_store_begin->first << "_" << component << "> ";
-                  }
-                         else
-                               to_gnuplot << "<" << data_store_begin->first << "_" << 0 << "> ";
+                unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components
+                std::vector <std::string> names = (component_names_map.find (data_store_begin->first))->second; 
+                if (names.size() > 0)
+                {
+                    for (unsigned int component = 0; component < names.size(); component++)
+                    {
+                        to_gnuplot << "<" << names[component] << "> ";
+                    }
+                }
+                else
+                {  
+                    for (unsigned int component = 0; component < n_comps; component++)
+                    {
+                        to_gnuplot << "<" << data_store_begin->first << "_" << component << "> ";
+                    }
+                }
             }
           to_gnuplot << "\n";
 
-                                          // write data stored for the point
+                       // write data stored for the point
           for (unsigned int key = 0; key < dataset_key.size (); key++)
             {
               std::map <std::string, std::vector <std::vector <double> > >::iterator
-               data_store_begin = data_store.begin ();
+                                              data_store_begin = data_store.begin ();
               to_gnuplot << dataset_key[key];
 
               for (unsigned int component = 0; component < n_indep; component++)
@@ -935,18 +1008,12 @@ void PointValueHistory<dim>
 
               for (; data_store_begin != data_store.end (); data_store_begin++)
                 {
-                                 bool scalar = (scalar_field.find (data_store_begin->first))->second;
-                                 if (scalar)
-                                 {
-                                       to_gnuplot << " " << (data_store_begin->second)[data_store_index][key];
-                                 }
-                                 else
-                                 {
-                                       for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
-                                       {
-                                         to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key];
-                                       }
-                                 }
+                  unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components
+                  
+                  for (unsigned int component = 0; component < n_comps; component++)
+                  {
+                      to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_comps + component][key];
+                  }
                 }
               to_gnuplot << "\n";
             }
@@ -962,11 +1029,11 @@ template <int dim>
 Vector<double> PointValueHistory<dim>
 ::mark_support_locations ()
 {
-                                  // a method to put a one at each point on
-                                  // the grid where a location is defined
-  AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
+                   // a method to put a one at each point on
+                   // the grid where a location is defined
   AssertThrow (!cleared, ExcInvalidState ());
-  AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+  AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
 
   Vector<double> dof_vector (dof_handler->n_dofs ());
 
@@ -986,7 +1053,7 @@ template <int dim>
 Vector<double> PointValueHistory<dim>
 ::mark_locations ()
 {
-       return mark_support_locations ();
+    return mark_support_locations ();
 }
 
 
@@ -994,9 +1061,9 @@ template <int dim>
 void PointValueHistory<dim>
 ::get_support_locations (std::vector <std::vector<Point <dim> > > & locations)
 {
-  AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
   AssertThrow (!cleared, ExcInvalidState ());
-  AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
+  AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
 
   std::vector <std::vector <Point <dim> > > actual_points;
   typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
@@ -1013,16 +1080,16 @@ template <int dim>
 void PointValueHistory<dim>
 ::get_points (std::vector <std::vector<Point <dim> > > & locations)
 {
-       get_support_locations (locations);
+    get_support_locations (locations);
 }
 
 
 template <int dim>
 void PointValueHistory<dim>
-::get_postprocessor_locations (std::vector<Point <dim> > & locations, const Quadrature<dim> & quadrature)
+::get_postprocessor_locations (const Quadrature<dim> & quadrature, std::vector<Point <dim> > & locations)
 {
   Assert (!cleared, ExcInvalidState ());
-  Assert (n_dofs != 0, ExcDoFHandlerRequired ());
+  AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
   
   locations = std::vector<Point <dim> > ();
 
@@ -1030,29 +1097,30 @@ void PointValueHistory<dim>
   unsigned int n_quadrature_points = quadrature.size();
   std::vector<Point<dim> > evaluation_points;
 
-// Loop over points and find correct cell
+                           // Loop over points and find correct cell
   typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
   for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
   {
-       // we now have a point to query, need to know what cell it is in
-       Point <dim> requested_location = point->requested_location;
-       typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1<dim>(), *dof_handler, requested_location).first;
-       fe_values.reinit (cell);
-       
-       evaluation_points = fe_values.get_quadrature_points();
-       double distance = cell->diameter ();
-       unsigned int selected_point = 0;
-       
-       for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
-       {
-         if (requested_location.distance (evaluation_points[q_point]) < distance)
-         {
-               selected_point = q_point;
-               distance = requested_location.distance (evaluation_points[q_point]);
-         }
-       }
-
-       locations.push_back (evaluation_points[selected_point]);
+                           // we now have a point to query, 
+                           // need to know what cell it is in
+    Point <dim> requested_location = point->requested_location;
+    typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1<dim>(), *dof_handler, requested_location).first;
+    fe_values.reinit (cell);
+    
+    evaluation_points = fe_values.get_quadrature_points();
+    double distance = cell->diameter ();
+    unsigned int selected_point = 0;
+    
+    for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
+    {
+      if (requested_location.distance (evaluation_points[q_point]) < distance)
+      {
+        selected_point = q_point;
+        distance = requested_location.distance (evaluation_points[q_point]);
+      }
+    }
+
+    locations.push_back (evaluation_points[selected_point]);
   }
 }
 
@@ -1064,6 +1132,8 @@ void PointValueHistory<dim>
   out << "***PointValueHistory status output***\n\n";
   out << "Closed: " << closed << "\n";
   out << "Cleared: " << cleared << "\n";
+  out << "Triangulation_changed: " << triangulation_changed << "\n";
+  out << "Have_dof_handler: " << have_dof_handler << "\n";
   out << "Geometric Data" << "\n";
 
   typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
@@ -1074,28 +1144,37 @@ void PointValueHistory<dim>
   else
     {
       if (!cleared)
-       {
-         for (; point != point_geometry_data.end (); point++)
-           {
-                 out << "# Requested location: " << point->requested_location << "\n";
-                 out << "# DoF_index : Support location (for each component)\n";
-             for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
-               {
-                 out << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n";
-               }
-             out << "\n";
-           }
-       }
+    {
+      for (; point != point_geometry_data.end (); point++)
+        {
+          out << "# Requested location: " << point->requested_location << "\n";
+          out << "# DoF_index : Support location (for each component)\n";
+          for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+        {
+          out << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n";
+        }
+          out << "\n";
+        }
+    }
       else
-       {
+    {
           out << "#Cannot access DoF_indices once cleared\n";
-       }
+    }
     }
   out << "\n";
 
   if (independent_values.size () != 0)
     {
       out << "Independent value(s): " << independent_values.size () << " : " << independent_values[0].size () << "\n";
+      if (indep_names.size() > 0)
+      {
+          out << "Names: ";
+          for (unsigned int name = 0; name < indep_names.size(); name++)
+          {
+              out << "<" << indep_names [name] << "> ";
+          }
+          out << "\n";
+      }
     }
   else
     {
@@ -1106,15 +1185,35 @@ void PointValueHistory<dim>
     data_store_begin = data_store.begin ();
   for (; data_store_begin != data_store.end (); data_store_begin++)
     {
+                  // Find field mnemonic 
+      std::string vector_name = data_store_begin->first;
+      typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
+      Assert (components != field_components.end(), ExcMessage("vector_name not in class"));
+      typename std::map <std::string, std::vector <std::string> >::iterator component_names = component_names_map.find(vector_name);
+      Assert (component_names != component_names_map.end(), ExcMessage("vector_name not in class"));
+      
       if (data_store_begin->second.size () != 0)
         {
-          out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << (data_store_begin->second)[0].size () << "\n";
+          out << data_store_begin->first << ": " << data_store_begin->second.size () << " (";
+          out << components->second.size() << ") : ";
+          out << (data_store_begin->second)[0].size () << "\n";
         }
       else
         {
-          out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << "No points added" << "\n";
+          out << data_store_begin->first << ": " << data_store_begin->second.size () << " (";
+          out << components->second.size() << ") : "; 
+          out << "No points added" << "\n";
         }
-    }
+                  // add names, if available
+      if (component_names->second.size() > 0)
+      {
+          for (unsigned int name = 0; name < component_names->second.size(); name++)
+          {
+            out << "<" << component_names->second[name] << "> ";
+          }
+          out << "\n";
+      }
+    } 
   out << "\n";
   out << "***end of status output***\n\n";
 }
@@ -1125,8 +1224,8 @@ template <int dim>
 bool PointValueHistory<dim>
 ::deep_check (bool strict)
 {
-                                  // test ways that it can fail, if control
-                                  // reaches last statement return true
+                   // test ways that it can fail, if control
+                   // reaches last statement return true
   if (strict)
     {
       if (n_indep != 0)
@@ -1137,19 +1236,19 @@ bool PointValueHistory<dim>
             }
         }
       std::map <std::string, std::vector <std::vector <double> > >::iterator
-       data_store_begin = data_store.begin ();
-      if (n_dofs != 0)
+    data_store_begin = data_store.begin ();
+      if (have_dof_handler)
         {
           for (; data_store_begin != data_store.end (); data_store_begin++)
             {
               if ((data_store_begin->second)[0].size () != dataset_key.size ())
                 return false;
-                                              // this loop only tests one
-                                              // member for each name,
-                                              // i.e. checks the user it will
-                                              // not catch internal errors
-                                              // which do not update all
-                                              // fields for a name.
+                           // this loop only tests one
+                           // member for each name,
+                           // i.e. checks the user it will
+                           // not catch internal errors
+                           // which do not update all
+                           // fields for a name.
             }
         }
       return true;
@@ -1162,19 +1261,19 @@ bool PointValueHistory<dim>
         }
     }
 
-  if (n_dofs != 0)
+  if (have_dof_handler)
     {
       std::map <std::string, std::vector <std::vector <double> > >::iterator
-       data_store_begin = data_store.begin ();
+    data_store_begin = data_store.begin ();
       for (; data_store_begin != data_store.end (); data_store_begin++)
         {
           if (std::abs ((int) (data_store_begin->second)[0].size () - (int) dataset_key.size ()) >= 2)
             return false;
-                                          // this loop only tests one member
-                                          // for each name, i.e. checks the
-                                          // user it will not catch internal
-                                          // errors which do not update all
-                                          // fields for a name.
+                       // this loop only tests one member
+                       // for each name, i.e. checks the
+                       // user it will not catch internal
+                       // errors which do not update all
+                       // fields for a name.
         }
     }
   return true;
@@ -1182,6 +1281,14 @@ bool PointValueHistory<dim>
 
 
 
+template <int dim>
+void PointValueHistory<dim>
+::tria_change_listener ()
+{
+  triangulation_changed = true;
+}
+
+
 // explicit instantiations
 #include "point_value_history.inst"
 
index 571307394b44005b489e1fdf44cea68ab95f69a2..e175ccbb234a3184765efbe17bcd00d761a16b1d 100644 (file)
@@ -23,8 +23,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS)
   template
   void PointValueHistory<deal_II_dimension>::evaluate_field 
   (const std::string &,
-   const VEC &,
-   const unsigned int);
+   const VEC &);
 }
 
 
@@ -33,15 +32,25 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS)
   template
   void PointValueHistory<deal_II_dimension>::evaluate_field_at_requested_location 
   (const std::string &,
+   const VEC &);
+}
+
+for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS)
+{  
+  template
+  void PointValueHistory<deal_II_dimension>::evaluate_field 
+  (const std::vector <std::string> &,
    const VEC &,
-   const unsigned int);
+   const DataPostprocessor<deal_II_dimension> &,
+   const Quadrature<deal_II_dimension> &);
 }
 
 for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS)
 {  
   template
   void PointValueHistory<deal_II_dimension>::evaluate_field 
-  (const VEC &,
+  (const std::string &,
+   const VEC &,
    const DataPostprocessor<deal_II_dimension> &,
    const Quadrature<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.