]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Michael Rapson: Major update of the PointValueHistory class.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 1 Jun 2012 19:07:07 +0000 (19:07 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 1 Jun 2012 19:07:07 +0000 (19:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@25595 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 9ab521fe121ef22597c806d18773e3070ecaf97f..5d3b023cfa8571fd434153f6dd349d24c72dc786 100644 (file)
@@ -26,6 +26,8 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/base/utilities.h>
+#include <numerics/data_postprocessor.h>
+#include <grid/grid_tools.h>
 
 #include <vector>
 #include <iostream>
@@ -59,9 +61,10 @@ namespace internal
     class PointGeometryData
     {
       public:
-        PointGeometryData(const std::vector <Point <dim> > &new_locations,
+        PointGeometryData(const Point <dim> &new_requested_location, const std::vector <Point <dim> > &new_locations,
                           const std::vector <int> &new_sol_indices);
-        std::vector <Point <dim> > locations;
+        Point <dim> requested_location;
+        std::vector <Point <dim> > support_point_locations;
         std::vector <int> solution_indices;
     };
   }
@@ -74,9 +77,10 @@ namespace internal
  * iterative process) graphs of solution values at specific points on the
  * mesh. The user specifies the points which the solution should be monitored
  * at ahead of time, as well as giving each solution vector that they want to
- * record a mnemonic name. Then for each step the user 'pushes back' the data
- * available from that time step and the class extracts data for the requested
- * points to store it. Finally once the computation is finished, the user can
+ * record a mnemonic name. Then, for each step the user calls one of the three
+ * available "evaluate field" methods to store the data from each time step,
+ * and the class extracts data for the requested
+ * points to store it. Finally, once the computation is finished, the user can
  * request Gnuplot output files to be generated.
  *
  * The user can store extra variables which do not relate to mesh location by
@@ -87,16 +91,44 @@ namespace internal
  * taken to solve the step and solver steps before convergence, or saving
  * norms calculated.
  *
- * Currently the code selects the nearest support points to a given point to
- * extract data from. This makes the code run at each time step very short,
+ * The three "evaluate field" methods each have different strengths and
+ * weaknesses making each suitable for different contexts:
+ * <ol>
+ * <li>Firstly, the @p evaluate_field version that does not take a @p DataPostprocessor object
+ * selects the nearest support points to a given point to
+ * extract data from. This makes the code that needs to be run at each time step very short,
  * since looping over the mesh to extract the needed dof_index can be done
- * just once at the start. However this does lead to problems when
- * FiniteElements which do not assign dofs to actual mesh locations are used
- * (i.e. FEs without support points). The Functions::FEFieldFunction class
- * allows users to evaluate a solution at any point in a domain and will work
- * even for FEs without support points. This functionality is not currently
- * offered through this class since it has more overhead per iteration than
- * the approach implemented, but this class can be extended to offer it.
+ * 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.)
+ *
+ * <li> Secondly, @p evaluate_field_at_requested_location calls @p
+ * VectorTools::point_value to compute values at the specific point requested.
+ * This method is valid for any FE that is supported by @p VectorTools::point_value.
+ * Specifically, this method can be called by codes using adaptive mesh refinement.
+ *
+ * <li>Finally, the class offers a function @p evaluate_field that takes a
+ * @p DataPostprocessor object. This method allows the deal.II
+ * data postprocessor to be used to compute new quantities from the
+ * solution on the fly. The values are located at the nearest quadrature
+ * point to the requested point. If the mesh is refined between calls, this
+ * point will change, so care must be taken when using this method in
+ * 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.
+ * </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.
  *
  * Following is a little code snippet that shows a common usage of this class:
  *
@@ -236,8 +268,8 @@ class PointValueHistory
                                       * point is found and its details stored
                                       * for use when @p evaluate_field is
                                       * called. If more than one point is
-                                      * required rather call this method as it
-                                      * is more efficient than the add_points
+                                      * required, rather call this method as it
+                                      * is more efficient than the add_point
                                       * method since it minimizes iterations
                                       * over the mesh. The points are added to
                                       * the internal database in the order
@@ -253,26 +285,104 @@ class PointValueHistory
 
                                      /**
                                       * Put another mnemonic string (and hence
-                                      * @p VECTOR) into the class. This also
-                                      * adds extra entries for points that are
-                                      * already in the class, so @p
-                                      * add_field_name and @p add_points can
+                                      * @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.
+                                      * 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);
 
+                                     /**
+                                      * 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.
+                                      * 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);
+
                                      /**
                                       * Extract values at the stored points
                                       * from the VECTOR supplied and add them
-                                      * to the new dataset in vector_name. If
-                                      * a @p DoFHandler is used, this method
+                                      * 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
+                                      * one) evaluate_field method
+                                      * must be called for each dataset (time
+                                      * step, iteration, etc) for each
+                                      * vector_name, otherwise a @p
+                                      * ExcDataLostSync error can occur.
+                                      */
+    template <class VECTOR>
+    void evaluate_field(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0);
+
+
+                                     /**
+                                      * Compute values using a DataPostprocessor object
+                                      * with the 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 quadrature object supplied is used for all
+                                      * components of a vector field.
+                                      * 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
+                                      * vector_name, otherwise a @p
+                                      * ExcDataLostSync error can occur.
+                                      */
+    template <class VECTOR>
+    void evaluate_field(const VECTOR & solution, const DataPostprocessor<dim> & data_postprocessor, const Quadrature<dim> & quadrature);
+
+
+                                     /**
+                                      * Extract values at the points actually
+                                      * requested from the VECTOR supplied and
+                                      * add them to the new dataset in vector_name.
+                                      * Unlike the other evaluate_field methods
+                                      * this method does not care if the dof_handler
+                                      * 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
+                                      * a @p DoFHandler is used, one (and only
+                                      * one) evaluate_field method
                                       * must be called for each dataset (time
                                       * step, iteration, etc) for each
                                       * vector_name, otherwise a @p
                                       * ExcDataLostSync error can occur.
                                       */
     template <class VECTOR>
-    void evaluate_field(const std::string &vector_name, const VECTOR & solution);
+    void evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0);
 
                                      /**
                                       * Add the key for the current dataset to
@@ -313,8 +423,23 @@ class PointValueHistory
                                       * key and independent data. The file
                                       * name agrees with the order the points
                                       * were added to the class.
+                                      * The support point information is only
+                                      * meaningful if the dof_handler has not
+                                      * been changed. Therefore, if adaptive
+                                      * mesh refinement has been used the
+                                      * support point data should not be used.
+                                      * The optional parameter
+                                      * postprocessor_locations is used to
+                                      * add the postprocessor locations to the
+                                      * output files. If this is desired, the
+                                      * data should be obtained from a call to
+                                      * get_postprocessor_locations while the
+                                      * dof_handler is usable. The default
+                                      * parameter is an empty vector of strings,
+                                      * and will suppress postprocessor
+                                      * locations output.
                                       */
-    void write_gnuplot(const std::string &base_name);
+    void write_gnuplot(const std::string &base_name, const std::vector <Point <dim> > postprocessor_locations = std::vector <Point <dim> > ());
 
 
                                      /**
@@ -336,15 +461,28 @@ class PointValueHistory
                                       * data_out.attach_dof_handler(dof_handler);
                                       *
                                       * // Call the mark_locations method to get the vector with indices flagged
-                                      * Vector<double> node_locations = node_monitor.mark_locations();
+                                      * Vector<double> support_point_locations = node_monitor.mark_locations();
                                       *
                                       * // Add the vector to the data_out object and write out a file in the usual way
-                                      * data_out.add_data_vector(node_locations, "Monitor_Locations");
+                                      * data_out.add_data_vector(support_point_locations, "Monitor_Locations");
                                       * data_out.build_patches(2);
                                       * std::ofstream output("locations.gpl");
                                       * data_out.write_gnuplot(output);
                                       * @endcode
                                       */
+    Vector<double> mark_support_locations();
+
+
+                                     /**
+                                      * Depreciated:
+                                      *
+                                      * This function only exists for backward
+                                      * compatibility as this is the interface
+                                      * provided by previous versions of the library.
+                                      * The function mark_support_locations replaces
+                                      * it and reflects the fact that the locations
+                                      * marked are actually the support points.
+                                      */
     Vector<double> mark_locations();
 
 
@@ -359,8 +497,37 @@ class PointValueHistory
                                       * the correct number of points by the
                                       * method.
                                       */
+    void get_support_locations (std::vector <std::vector<Point <dim> > > & locations);
+
+                                     /**
+                                      * Depreciated:
+                                      *
+                                      * This function only exists for backward
+                                      * compatibility as this is the interface
+                                      * provided by previous versions of the library.
+                                      * The function get_support_locations replaces
+                                      * it and reflects the fact that the points
+                                      * returned are actually the support points.
+                                      */
     void get_points (std::vector <std::vector<Point <dim> > > & locations);
 
+                                     /**
+                                      * Stores the actual location of the
+                                      * points used by the data_postprocessor.
+                                      * This can be used to compare with the
+                                      * points requested, for example by using
+                                      * the @p Point<dim>::distance function.
+                                      * Unlike the support_locations, these
+                                      * locations are computed every time the
+                                      * evaluate_field method is called with a
+                                      * postprocessor. This method uses the
+                                      * same algorithm so can will find the
+                                      * same points.
+                                      * For convenience, location is resized
+                                      * to the correct number of points by the
+                                      * method.
+                                      */
+    void get_postprocessor_locations (std::vector<Point <dim> > & locations, const Quadrature<dim> & quadrature);
 
                                      /**
                                       * Once datasets have been added to the
@@ -490,11 +657,24 @@ class PointValueHistory
 
                                      /**
                                       * Saves data for each mnemonic entry.
-                                      * data_store: mnemonic -> [component]
-                                      * [key]
+                                      * data_store: mnemonic ->
+                                      * [point_0_components point_1_components
+                                      * ... point_n-1_components][key]
+                                      * This format facilitates scalar mnemonics
+                                      * in a vector space, because scalar mnemonics
+                                      * will only have one component per point.
+                                      * Vector components are strictly
+                                      * 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).
+                                      */
+    std::map <std::string, bool> scalar_field;
+
                                      /**
                                       * Saves the location and other mesh info
                                       * about support points.
@@ -508,43 +688,42 @@ class PointValueHistory
                                       */
     std::pair<std::string, std::vector <std::vector <double> > > pair_data; // could possibly be removed
 
-                                    /**
-                                     * Used to enforce @p closed state for some
-                                     * methods.
-                                     */
+                                     /**
+                                      * Used to enforce @p closed state for some
+                                      * methods.
+                                      */
     bool closed;
 
-                                    /**
-                                     * Used to enforce @p !cleared state for
-                                     * some methods.
-                                     */
+                                     /**
+                                      * Used to enforce @p !cleared state for
+                                      * some methods.
+                                      */
     bool cleared;
 
 
-                                    /**
-                                     * A smart pointer to the dof_handler
-                                     * supplied to the constructor. This can be
-                                     * released by calling @p clear().
-                                     */
+                                     /**
+                                      * A smart pointer to the dof_handler
+                                      * supplied to the constructor. This can be
+                                      * released by calling @p clear().
+                                      */
     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
+                                      *
+                                     * 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.
                                       */
     unsigned int n_dofs;
-                                    /**
-                                     * Stores the number of independent
-                                     * variables requested.
-                                     */
+                                     /**
+                                      * Stores the number of independent
+                                      * variables requested.
+                                      */
     unsigned int n_indep;
-
 };
 
 
index a1ea25af9703657f28b30ddb0fa265133830f2f3..2cfa0ea914ed8cf6354da78b9cb40ecf2d0a2403 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2009, 2010, 2012 by Michael Rapson and the deal.II authors
+//    Copyright (C) 2009, 2010 by Michael Rapson and the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -20,6 +20,8 @@
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/trilinos_block_vector.h>
 
+#include <deal.II/numerics/vectors.h>
+
 #include <deal.II/numerics/point_value_history.h>
 
 
@@ -33,10 +35,12 @@ namespace internal
 /// Only a constructor needed for this class (a struct really)
     template <int dim>
     PointGeometryData<dim>
-    ::PointGeometryData (const std::vector <Point <dim> > &new_locations,
-                         const std::vector <int> &new_sol_indices)
+    ::PointGeometryData (const Point <dim> &new_requested_location, 
+                        const std::vector <Point <dim> > &new_locations,
+                        const std::vector <int> &new_sol_indices)
     {
-      locations = new_locations;
+         requested_location = new_requested_location;
+      support_point_locations = new_locations;
       solution_indices = new_sol_indices;
     }
   }
@@ -47,15 +51,15 @@ namespace internal
 template <int dim>
 PointValueHistory<dim>
 ::PointValueHistory (const unsigned int n_independent_variables) :
-                n_dofs (0),
-                n_indep (n_independent_variables)
+               n_dofs (0),
+               n_indep (n_independent_variables)
 {
   closed = false;
   cleared = false;
-                                   // make a vector for keys
+                                  // 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));
 }
@@ -64,17 +68,17 @@ PointValueHistory<dim>
 
 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_dofs (dof_handler.n_dofs ()),
+               n_indep (n_independent_variables)
 {
   closed = false;
   cleared = false;
-                                   // make a vector to store keys
+                                  // 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));
 }
@@ -87,6 +91,7 @@ PointValueHistory<dim>::PointValueHistory (const PointValueHistory & point_value
   dataset_key = point_value_history.dataset_key;
   independent_values = point_value_history.independent_values;
   data_store = point_value_history.data_store;
+  scalar_field = point_value_history.scalar_field;
   point_geometry_data = point_value_history.point_geometry_data;
   pair_data = point_value_history.pair_data;
   closed = point_value_history.closed;
@@ -107,6 +112,7 @@ PointValueHistory<dim>::operator= (const PointValueHistory & point_value_history
   dataset_key = point_value_history.dataset_key;
   independent_values = point_value_history.independent_values;
   data_store = point_value_history.data_store;
+  scalar_field = point_value_history.scalar_field;
   point_geometry_data = point_value_history.point_geometry_data;
   pair_data = point_value_history.pair_data;
   closed = point_value_history.closed;
@@ -133,49 +139,49 @@ 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 ()));
+              ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
 
-                                   // 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);
@@ -183,32 +189,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;
@@ -221,44 +227,45 @@ 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>
-    new_point_geometry_data (current_points, new_solution_indices);
+    new_point_geometry_data (location, current_points, new_solution_indices);
   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 ();
   for (; data_store_begin != data_store.end (); data_store_begin++)
     {
-                                       // add an extra row to each vector
-                                       // entry
+                                      // add an extra row to each vector
+                                      // entry
+         bool scalar = (scalar_field.find (data_store_begin->first))->second;
       for (unsigned int component = 0;
-           component < dof_handler->get_fe ().n_components (); component++)
+          component < (scalar ? 1 : dof_handler->get_fe ().n_components ()); component++)
         {
           data_store_begin->second.push_back (std::vector<double> (0));
         }
@@ -271,46 +278,46 @@ 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
+                                  // 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);
@@ -318,8 +325,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;
@@ -327,8 +334,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);
@@ -341,7 +348,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;
@@ -361,17 +368,18 @@ void PointValueHistory<dim>
           new_solution_indices.push_back (local_dof_indices[current_fe_index[point][component]]);
         }
 
-      internal::PointValueHistory::PointGeometryData<dim> new_point_geometry_data (current_points[point], new_solution_indices);
+      internal::PointValueHistory::PointGeometryData<dim> new_point_geometry_data (locations[point], current_points[point], new_solution_indices);
 
       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
-          for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+                                          // 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++)
             {
               data_store_begin->second.push_back (std::vector<double> (0));
             }
@@ -388,22 +396,53 @@ template <int dim>
 void PointValueHistory<dim>
 ::add_field_name (const std::string &vector_name)
 {
-                                   // 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 ()));
 
 
-                                   // make and add a new vector
-                                   // point_geometry_data.size() long
+                                  // make and add a new vector
+                                  // point_geometry_data.size() long
   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
   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)
+{
+                                  // 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 ()));
+
+
+                                  // 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
 }
 
 
@@ -440,10 +479,10 @@ void PointValueHistory<dim>
 template <int dim>
 template <class VECTOR>
 void PointValueHistory<dim>
-::evaluate_field (const std::string &vector_name, const VECTOR & solution)
+::evaluate_field (const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component)
 {
-                                   // 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 ());
@@ -452,28 +491,297 @@ void PointValueHistory<dim>
     {
       Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
     }
-  typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
-  for (unsigned int data_store_index = 0; node != point_geometry_data.end (); node++, data_store_index++)
+                                  // 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"));
+  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++)
     {
-                                       // step through each node, 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 = node->solution_indices[component];
-          (data_store[vector_name])[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (solution (solution_index));
-        }
+         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));
+         }
     }
 }
 
 
 
+
+
+template <int dim>
+template <class VECTOR>
+void PointValueHistory<dim>
+::evaluate_field(const VECTOR & solution, const DataPostprocessor< dim> & data_postprocessor, const Quadrature<dim> & quadrature)
+{
+                                  // must be closed to add data to internal
+                                  // members.
+  Assert (closed, ExcInvalidState ());
+  Assert (!cleared, ExcInvalidState ());
+  Assert (n_dofs != 0, 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
+  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."));
+  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 ();
+
+// 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 ++;
+                       }
+               }
+    } // end of loop over points
+}
+
+
+
+
+
+
+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)
+{
+                                  // must be closed to add data to internal
+                                  // members.
+  Assert (closed, ExcInvalidState ());
+  Assert (!cleared, ExcInvalidState ());
+  Assert (n_dofs != 0, 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.
+  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"));
+  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);
+         
+         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));
+         }
+    }
+}
+
+
 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 ());
@@ -487,8 +795,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));
@@ -503,13 +811,13 @@ void PointValueHistory<dim>
 
 template <int dim>
 void PointValueHistory<dim>
-::write_gnuplot (const std::string &base_name)
+::write_gnuplot (const std::string &base_name, const std::vector <Point <dim> > postprocessor_locations)
 {
   AssertThrow (closed, ExcInvalidState ());
   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";
@@ -517,7 +825,7 @@ 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++)
@@ -526,7 +834,7 @@ void PointValueHistory<dim>
         }
       to_gnuplot << "\n";
 
-                                       // write general data stored
+                                      // write general data stored
       for (unsigned int key = 0; key < dataset_key.size (); key++)
         {
           to_gnuplot << dataset_key[key];
@@ -543,39 +851,57 @@ void PointValueHistory<dim>
 
 
 
-                                   // write points to a file
+                                  // write points to a file
   if (n_dofs != 0)
     {
       AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
-      AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+         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 ();
 
-      typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
-      for (unsigned int data_store_index = 0; node != point_geometry_data.end (); node++, data_store_index++)
+      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 << "# DoF_index : Location (for each component)\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++)
             {
-              to_gnuplot << "# " << node->solution_indices[component] << " : " << node->locations [component] << "\n";
+              to_gnuplot << "# " << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n";
             }
-          to_gnuplot << "\n";
+                 if (postprocessor_locations.size() != 0)
+                 {
+                         to_gnuplot << "# Postprocessor location: " << postprocessor_locations[data_store_index] << "\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++)
@@ -585,18 +911,21 @@ void PointValueHistory<dim>
 
           for (; data_store_begin != data_store.end (); data_store_begin++)
             {
-              for (unsigned int component = 0; component < n_components; component++)
-                {
-                  to_gnuplot << "<" << data_store_begin->first << "_" << component << "> ";
-                }
+                         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 << "> ";
             }
           to_gnuplot << "\n";
 
-                                           // write data stored for the node
+                                          // 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++)
@@ -606,10 +935,18 @@ void PointValueHistory<dim>
 
               for (; data_store_begin != data_store.end (); data_store_begin++)
                 {
-                  for (unsigned int component = 0; component < n_components; component++)
-                    {
-                      to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key];
-                    }
+                                 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];
+                                       }
+                                 }
                 }
               to_gnuplot << "\n";
             }
@@ -623,48 +960,102 @@ void PointValueHistory<dim>
 
 template <int dim>
 Vector<double> PointValueHistory<dim>
-::mark_locations ()
+::mark_support_locations ()
 {
-                                   // a method to put a one at each point on
-                                   // the grid where a location is defined
+                                  // a method to put a one at each point on
+                                  // the grid where a location is defined
   AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
   AssertThrow (!cleared, ExcInvalidState ());
   AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
 
   Vector<double> dof_vector (dof_handler->n_dofs ());
 
-  typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
-  for (; node != point_geometry_data.end (); node++)
+  typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+  for (; point != point_geometry_data.end (); point++)
     {
       for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
         {
-          dof_vector (node->solution_indices[component]) = 1;
+          dof_vector (point->solution_indices[component]) = 1;
         }
     }
   return dof_vector;
 }
 
 
+template <int dim>
+Vector<double> PointValueHistory<dim>
+::mark_locations ()
+{
+       return mark_support_locations ();
+}
+
 
 template <int dim>
 void PointValueHistory<dim>
-::get_points (std::vector <std::vector<Point <dim> > > & locations)
+::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 ()));
 
   std::vector <std::vector <Point <dim> > > actual_points;
-  typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
+  typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
 
-  for (; node != point_geometry_data.end (); node++)
+  for (; point != point_geometry_data.end (); point++)
     {
-      actual_points.push_back (node->locations);
+      actual_points.push_back (point->support_point_locations);
     }
   locations = actual_points;
 }
 
 
+template <int dim>
+void PointValueHistory<dim>
+::get_points (std::vector <std::vector<Point <dim> > > & locations)
+{
+       get_support_locations (locations);
+}
+
+
+template <int dim>
+void PointValueHistory<dim>
+::get_postprocessor_locations (std::vector<Point <dim> > & locations, const Quadrature<dim> & quadrature)
+{
+  Assert (!cleared, ExcInvalidState ());
+  Assert (n_dofs != 0, ExcDoFHandlerRequired ());
+  
+  locations = std::vector<Point <dim> > ();
+
+  FEValues<dim> fe_values (dof_handler->get_fe (), quadrature, update_quadrature_points);
+  unsigned int n_quadrature_points = quadrature.size();
+  std::vector<Point<dim> > evaluation_points;
+
+// 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]);
+  }
+}
+
 
 template <int dim>
 void PointValueHistory<dim>
@@ -675,29 +1066,30 @@ void PointValueHistory<dim>
   out << "Cleared: " << cleared << "\n";
   out << "Geometric Data" << "\n";
 
-  typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
-  if (node == point_geometry_data.end ())
+  typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+  if (point == point_geometry_data.end ())
     {
-      out << "No nodes stored currently\n";
+      out << "No points stored currently\n";
     }
   else
     {
       if (!cleared)
-        {
-          out << "# DoF_index : Location (for each component)\n";
-          for (; node != point_geometry_data.end (); node++)
-            {
-              for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
-                {
-                  out << node->solution_indices[component] << " : " << node->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";
 
@@ -733,8 +1125,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)
@@ -745,19 +1137,19 @@ bool PointValueHistory<dim>
             }
         }
       std::map <std::string, std::vector <std::vector <double> > >::iterator
-        data_store_begin = data_store.begin ();
+       data_store_begin = data_store.begin ();
       if (n_dofs != 0)
         {
           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;
@@ -773,16 +1165,16 @@ bool PointValueHistory<dim>
   if (n_dofs != 0)
     {
       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;
index 9d3fdbeaf1651d66f5e64280c6ba93adb277ecf9..571307394b44005b489e1fdf44cea68ab95f69a2 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2009, 2012 by Michael Rapson and the deal.II authors
+//    Copyright (C) 2009 by Michael Rapson and the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -19,9 +19,29 @@ for (deal_II_dimension : DIMENSIONS)
 
 
 for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS)
-{
+{  
   template
-  void PointValueHistory<deal_II_dimension>::evaluate_field
+  void PointValueHistory<deal_II_dimension>::evaluate_field 
   (const std::string &,
-   const VEC &);
+   const VEC &,
+   const unsigned int);
+}
+
+
+for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS)
+{  
+  template
+  void PointValueHistory<deal_II_dimension>::evaluate_field_at_requested_location 
+  (const std::string &,
+   const VEC &,
+   const unsigned int);
+}
+
+for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS)
+{  
+  template
+  void PointValueHistory<deal_II_dimension>::evaluate_field 
+  (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.