]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Adding point_value_history class to numerics and compute_projection_from_face_quadrat...
authorrapson <rapson@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 Jul 2009 15:43:50 +0000 (15:43 +0000)
committerrapson <rapson@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 8 Jul 2009 15:43:50 +0000 (15:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@19042 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_tools.h
deal.II/deal.II/include/numerics/point_value_history.h [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_tools.cc
deal.II/deal.II/source/numerics/point_value_history.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/point_value_history.inst.in [new file with mode: 0644]

index 9eebc513798f6e711e17fb8087be621f1203f4d7..6fae20c4f43bb7f253dcf802727af84e4db42273 100644 (file)
@@ -644,6 +644,31 @@ class FETools
     compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
                                                        const Quadrature<dim>    &quadrature,
                                                        FullMatrix<double>       &I_q);
+
+
+
+
+                                    /**
+                                     * This method implements the 
+                                     * FETools::compute_projection_from_quadrature_points_matrix
+                                     * method for faces of a mesh. 
+                                     * The matrix that it returns, X, is face specific
+                                     * and its size is fe.dofs_per_cell by
+                                     * rhs_quadrature.size().
+                                     * The dimension, dim must be larger than 1 for this class,
+                                     * since Quadrature<dim-1> objects are required. See the 
+                                     * documentation on the Quadrature class for more information.
+                                     */
+    template <int dim, int spacedim>
+    static
+    void
+    compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
+                                               const Quadrature<dim-1>    &lhs_quadrature,
+                                               const Quadrature<dim-1>    &rhs_quadrature,
+                                               const typename DoFHandler<dim, spacedim>::active_cell_iterator & cell,
+                                               unsigned int face,
+                                               FullMatrix<double>       &X);
+
     
 
                                     //@}
@@ -1330,6 +1355,17 @@ class FETools
                                      */
     DeclException1(ExcLeastSquaresError, double,
                   << "Least squares fit leaves a gap of " << arg1);
+
+                                    /**
+                                     * Exception thrown if one variable
+                                     * may not be greater than another.
+                                     *
+                                     * @ingroup Exceptions
+                                     */
+    DeclException2 (ExcNotGreaterThan,
+                  int, int,
+                  << arg1 << " must be greater than " << arg2);
+
   private:
                                     /**
                                      * Return a finite element that
diff --git a/deal.II/deal.II/include/numerics/point_value_history.h b/deal.II/deal.II/include/numerics/point_value_history.h
new file mode 100644 (file)
index 0000000..45f392c
--- /dev/null
@@ -0,0 +1,365 @@
+//---------------------------------------------------------------------------
+//
+//    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
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __dealii__point_value_history_h
+#define __dealii__point_value_history_h
+
+#include <dofs/dof_handler.h>
+#include <base/point.h>
+#include <lac/vector.h>
+
+#include <dofs/dof_accessor.h>
+#include <base/exceptions.h>
+#include <base/quadrature_lib.h>
+#include <fe/fe_q.h>
+#include <fe/mapping.h>
+#include <fe/mapping_q1.h>
+#include <fe/fe_values.h>
+#include <base/smartpointer.h>
+#include <base/utilities.h>
+
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <map>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+    namespace PointValueHistory
+    {
+        template <int dim> class PointGeometryData;
+    }
+}
+
+
+
+namespace internal
+{
+    namespace PointValueHistory
+    {    
+        /*!
+         *  A class that stores the data needed to reference the support points closest to one
+         *  requested point.
+         */
+        template <int dim>
+        class PointGeometryData
+        {
+        public:
+            PointGeometryData(std::vector <Point <dim> > new_locations, std::vector <int> new_sol_indices);
+            std::vector <Point <dim> > locations;
+            std::vector <int> solution_indices;
+        };
+    }
+}
+
+
+
+/*!
+ * PointValueHistory tackles the overhead of plotting time (or any other 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
+ * request Gnuplot output files to be generated.
+ *
+ * The user can store extra variables which do not relate to mesh location by the parameter by specifying
+ * n_independent_variables. The class then expects a std::vector of size n_independent_variables to be added
+ * during each step using the method @p push_back_independent. This may be used for example for recording 
+ * external input, logging solver performance data such as time 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, 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, (FEs without support points). The FEFieldFunction class
+ * allows users to evaluate a solution at any point in dof handler 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.
+ *
+ * Code snippets showing common usage of the class flows:
+ *
+ * @code
+ * #include <numerics/point_value_history.h>
+ * //....
+ *
+ * //... code to setup Triangulation, perform any refinement necessary
+ * // and setup DoFHandler, sizing solution Vectors etc
+ *
+ * // call the constructor
+ * unsigned int n_inputs = 1; // just one independent value, which happens to be an input
+ * PointValueHistory<dim> node_monitor(dof_handler, n_inputs);
+ * // or if node_monitor_2 is defined in class
+ * node_monitor_2 = PointValueHistory<dim> (dof_handler, n_inputs);
+ *
+ * // setup fields and points required
+ * node_monitor.add_field_name("Solution");
+ * std::vector <Point <dim> > point_vector(2, Point <dim> ::Point());
+ * point_vector[0] = Point <dim> ::Point(0, 0);
+ * point_vector[1] = Point <dim> ::Point(0.25, 0);
+ * node_monitor.add_points(point_vector); // multiple points at once
+ * node_monitor.add_point(Point<dim>::Point(1, 0.2)); // add a single point
+ * node_monitor.close(); // close the class once the setup is complete
+ * node_monitor.status(std::cout); // print out status to check if desired
+ *
+ * // ... more code ...
+ *
+ * // ... in an iterative loop ...
+ * // double time, vector <double> with size 1 input_value,
+ * // and Vector <double> solution calculated in the loop
+ * node_monitor.start_new_dataset(time);
+ * node_monitor.push_back_independent(input_value);
+ * node_monitor.evaluate_field("Solution", solution);
+ *
+ * // ... end of iterative loop ...
+ *
+ * node_monitor.write_gnuplot("node"); // write out data files
+ * node_monitor_2.clear (); //clear dof_handler (only if required)
+ *
+ * @endcode
+ */
+template <int dim>
+class PointValueHistory
+{
+public:
+    /*!
+     * Provide a stripped down instance of the class which does not support adding points or mesh data.
+     * This may be used for example for recording external input or logging solver performance data.
+     */
+    PointValueHistory (unsigned int n_independent_variables = 0);
+    
+    /*!
+     * Constructor linking the class to a specific @p DoFHandler. Since 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.
+     * Triangulation 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.
+     *
+     * The user can store extra variables which do not relate to mesh location by specifying the 
+     * number required using n_independent_variables and making calls to @p push_back_independent as needed. 
+     * This may be used for example for recording external input or logging solver performance data.
+     */
+    PointValueHistory (const DoFHandler<dim> & dof_handler, unsigned int n_independent_variables = 0);
+    /*!
+     * Copy constructor explicitly provided, although default copy operator would be sufficient.
+     * This constructor can be safely called with a @p PointValueHistory object that contains data, but
+     * this could be expensive and should be avoided.
+     */
+    PointValueHistory (const PointValueHistory & point_value_history);
+
+    /*!
+     * Assignment operator explicitly provided, although default assignment operator would be sufficient.
+     * This assignment operator can be safely called once the class is closed and data added, but this
+     * is provided primarily to allow a @p PointValueHistory object declared in a class to be reinitialized
+     * later in the class. Using the assignment operator when the object contains data could be expensive.
+    */
+    PointValueHistory & operator=(const PointValueHistory & point_value_history);
+
+    /*!
+    Deconstructor, nothing actually required!
+    */
+    ~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 for use when @p evaluate_field is called. If more than one point is required
+     * rather used the @p add_points method since this minimizes iterations over the mesh.
+     */
+     void add_point(const Point <dim> & location);
+      
+    /*!
+     * 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 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 method since it minimizes 
+     * iterations over the mesh. The points are added to the internal database in the order they appear in the
+     * list and there is always a one to one correspondence between the requested point and the added point,
+     * even if a point is requested multiple times.
+     */
+    void add_points (const std::vector <Point <dim> > & locations);
+    
+
+      
+    /*!
+     * 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 be called in any order.
+     */
+    void add_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 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);
+    
+    /*!
+     * Add the key for the current dataset to the dataset. Although calling this method first is 
+     * sensible, the order in which this method, @p evaluate_field and @p push_back_independent is
+     * not important. It is however important that all the data for a give dataset is added to each
+     * dataset and that it is added before a new data set is started. This prevents a @p ExcDataLostSync.
+     */
+    void start_new_dataset(double key);
+    
+    /*!
+     * If independent values have been set up, this method stores these values. This should only be called
+     * once per dataset, and if independent values are used it must be called for every dataset. A @p ExcDataLostSync
+     * exception can be thrown if this method is not called.
+     */
+    void push_back_independent(const std::vector <double> &independent_values);
+    
+    
+    /*!
+     * Write out a series of .gpl files named base_name + "-00.gpl", base_name + "-01.gpl" etc. The
+     * data file gives info about where the support points selected and interpreting the data. 
+     * If @p n_indep != 0 an additional file base_name + "_indep.gpl" containing key and independent 
+     * data. The file name agrees with the order the points were added to the class.
+     */
+    void write_gnuplot(const std::string &base_name);
+    
+    
+    /*!
+     * Return a @p Vector with the indices of selected points flagged with a 1. This method is mainly
+     * for testing and verifying that the class is working correctly. By passing this vector to a 
+     * DataOut object, the user can verify that the positions returned by @p PointValueHistory< dim >::get_points 
+     * agree with the positions that @p DataOut interprets from the @p Vector returned. The code snippet below
+     * demonstrates how this could be done:
+     * @code
+           // Make a DataOut object and attach the dof_handler 
+           DataOut<dim> data_out;
+           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();
+     *
+           // 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.build_patches(2);
+           std::ofstream output("locations.gpl");
+           data_out.write_gnuplot(output); @endcode
+     */
+    
+    Vector<double> mark_locations();
+    
+    
+    /*!
+     * Stores the actual location of each support point selected by the @p add_point(s) method.
+     * This can be used to compare with the point requested, for example by using the 
+     * @p Point<dim>::distance function. For convenience, location is resized to the 
+     * correct number of points by the method.
+     */
+     void get_points (std::vector <std::vector<Point <dim> > > & locations);
+    
+    
+    /*!
+     * Once datasets have been added to the class, requests to add additional points
+     * will make the data interpretation unclear. The boolean @p closed defines a state of the 
+     * class and ensures this does not happen. Additional points or vectors can only be added 
+     * while the class is not closed, and the class must be closed before datasets can be added or 
+     * written to file. @p PointValueHistory::get_points and @p PointValueHistory::status do not require
+     * the class to be closed. If a method that requires a class to be open or close is called while in
+     * the wrong state a @p ExcInvalidState exception is thrown.
+     */
+    void close();
+
+
+    /*!
+     * Delete the lock this object has to the @p DoFHandler used the last time the class was created.
+     * This method should not normally need to be called, but can be useful to ensure that the
+     * @p DoFHandler is released before it goes out of scope if the @p PointValueHistory class might
+     * live longer than it. Once this method has been called, the majority of methods will throw a
+     * @p ExcInvalidState exception, so if used this method should be the last call to the class.
+     */
+    void clear();
+    
+    /*!
+     * Print useful debugging information about the class, include details about which support 
+     * points were selected for each point and sizes of the data stored.
+     */
+    void status(std::ostream &out);
+    
+    
+    /*!
+     * Check the internal data sizes to test for a loss of data sync. This is often used in 
+     * @p Assert statements with the @p ExcDataLostSync exception. If @p strict is @p false this method 
+     * returns @p true if all sizes are within 1 of each other (needed to allow data to be added), 
+     * with @p strict = @p true they must be exactly equal.
+     */
+    
+    bool deep_check (bool strict);
+   
+    /*!
+     * A call has been made to @p push_back_independent when no independent values were requested. 
+     */
+    DeclException0(ExcNoIndependent);   
+
+    /*!
+     * This error is thrown to indicate that the data sets appear to be out of sync.
+     * The class requires that the number of dataset keys is the same as the number of independent
+     * values sets and mesh linked value sets. The number of each of these is allowed to differ by one to 
+     * allow new values to be added with out restricting the order the user choses to do so.
+     * Special cases of no @p DoFHandler and no independent values should not trigger this error.     
+     */
+    DeclException0(ExcDataLostSync);    
+    
+    
+    /*!
+     * 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 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 that
+     * the internal dof indices stored are no longer meaningful. 
+     */
+    DeclException2(ExcDoFHandlerChanged,
+            int,
+            int,
+            << "Original n_dofs (" << arg1 << ") != current n_dofs (" << arg2 << ")");
+
+private:
+
+    std::vector <double> dataset_key; // commonly time, but possibly time step, iteration etc
+    std::vector <std::vector <double> > independent_values; // values that do not depend on grid location
+    std::map <std::string, std::vector <std::vector <double> > > data_store; // Saves data for each mnemonic entry.  data_store: mnemonic -> [component] [time_instance]
+    std::vector <internal::PointValueHistory::PointGeometryData <dim> > point_geometry_data; // Saves the location and other mesh info about support points
+    std::pair<std::string, std::vector <std::vector <double> > > pair_data; // could possibly be removed, just used as the value to be added to data_store.
+    bool closed;
+    bool cleared; 
+
+    SmartPointer<const DoFHandler<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.
+     */
+    
+    unsigned int n_dofs;
+    unsigned int n_indep;
+
+};
+
+
+DEAL_II_NAMESPACE_CLOSE
+#endif /* __dealii__point_value_history_h */
index 9e2c36dbd28bce862999bc012ab8878f68a20143..dd34c061ca04a6db49106b8f62c4a0e9df7aa242 100644 (file)
@@ -1958,6 +1958,68 @@ compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,space
 
 
 
+template <int dim, int spacedim>
+void
+FETools::
+compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
+                                               const Quadrature<dim-1>    &lhs_quadrature,
+                                               const Quadrature<dim-1>    &rhs_quadrature,
+                                               const typename DoFHandler<dim, spacedim>::active_cell_iterator & cell,
+                                               unsigned int face,
+                                               FullMatrix<double>       &X)
+{
+  Assert (fe.n_components() == 1, ExcNotImplemented());
+  Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree));
+
+
+
+                                   // build the matrices M and Q
+                                   // described in the documentation
+  FullMatrix<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
+  FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
+
+  {
+                               // need an FEFaceValues object to evaluate shape function
+                               // values on the specified face.
+    FEFaceValues <dim> fe_face_values (fe, lhs_quadrature, update_values); 
+    fe_face_values.reinit (cell, face); // setup shape_value on this face.
+
+    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+      for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+        for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
+          M(i,j) += fe_face_values.shape_value (i, q) *
+                    fe_face_values.shape_value (j, q) *
+                    lhs_quadrature.weight(q);
+    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+    {
+      M(i,i) = (M(i,i) == 0 ? 1 : M(i,i));
+    }
+  }
+
+  {
+    FEFaceValues <dim> fe_face_values (fe, rhs_quadrature, update_values); 
+    fe_face_values.reinit (cell, face); // setup shape_value on this face.
+
+    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+      for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
+        Q(i,q) += fe_face_values.shape_value (i, q) *
+                  rhs_quadrature.weight(q);
+  }
+                                   // then invert M
+  FullMatrix<double> M_inverse (fe.dofs_per_cell, fe.dofs_per_cell);
+  M_inverse.invert (M);
+
+                                   // finally compute the result
+  X.reinit (fe.dofs_per_cell, rhs_quadrature.size());
+  M_inverse.mmult (X, Q);
+
+  Assert (X.m() == fe.dofs_per_cell, ExcInternalError());
+  Assert (X.n() == rhs_quadrature.size(), ExcInternalError());
+}
+
+
+
+
 /*-------------- Explicit Instantiations -------------------------------*/
 
 
@@ -2412,6 +2474,18 @@ compute_interpolation_to_quadrature_points_matrix (const FiniteElement<deal_II_d
                                                    const Quadrature<deal_II_dimension>    &quadrature,
                                                    FullMatrix<double>       &I_q);
 
+#if deal_II_dimension != 1
+template
+void
+FETools::
+compute_projection_from_face_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
+                                               const Quadrature<deal_II_dimension-1>    &lhs_quadrature,
+                                               const Quadrature<deal_II_dimension-1>    &rhs_quadrature,
+                                               const DoFHandler<deal_II_dimension>::active_cell_iterator & cell,
+                                               unsigned int face,
+                                               FullMatrix<double>       &X);
+#endif
+
 
 
 /*----------------------------   fe_tools.cc     ---------------------------*/
diff --git a/deal.II/deal.II/source/numerics/point_value_history.cc b/deal.II/deal.II/source/numerics/point_value_history.cc
new file mode 100644 (file)
index 0000000..129bcd5
--- /dev/null
@@ -0,0 +1,708 @@
+//---------------------------------------------------------------------------
+//
+//    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
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <numerics/point_value_history.h>
+
+
+
+DEAL_II_NAMESPACE_OPEN
+
+#include "point_value_history.inst"
+
+/// Only a constructor needed for this class (a struct really)
+template <int dim>
+internal::PointValueHistory::PointGeometryData<dim>
+::PointGeometryData (std::vector <Point <dim> > new_locations, std::vector <int> new_sol_indices)
+{
+  locations = new_locations;
+  solution_indices = new_sol_indices;
+}
+
+
+
+template <int dim>
+PointValueHistory<dim>
+::PointValueHistory (unsigned int n_independent_variables) :
+dof_handler (0, "No DoFHandler"),
+n_dofs (0),
+n_indep (n_independent_variables)
+{
+  closed = false;
+  cleared = false;
+  // make a vector for keys
+  dataset_key = std::vector <double> (); // initialize the std::vector
+
+  // make a vector of independent values
+  independent_values = std::vector<std::vector <double> > (n_indep, std::vector <double> (0));
+}
+
+
+
+template <int dim>
+PointValueHistory<dim>
+::PointValueHistory (const DoFHandler<dim> & dof_handler, unsigned int n_independent_variables) :
+dof_handler (&dof_handler, typeid (*this).name ()),
+n_dofs (dof_handler.n_dofs ()),
+n_indep (n_independent_variables)
+{
+  closed = false;
+  cleared = false;
+  // make a vector to store keys
+  dataset_key = std::vector <double> (); // initialize the std::vector
+
+  // make a vector for the independent values
+  independent_values = std::vector<std::vector <double> > (n_indep, std::vector <double> (0));
+}
+
+
+
+template <int dim>
+PointValueHistory<dim>
+::PointValueHistory (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;
+  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;
+  n_indep = point_value_history.n_indep;
+}
+
+
+
+template <int dim>
+        PointValueHistory<dim> & 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;
+  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;
+  n_indep = point_value_history.n_indep;
+  return * this;
+}
+
+
+
+template <int dim>
+PointValueHistory<dim>
+::~PointValueHistory ()
+{
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::add_point (const Point <dim> & location)
+{
+  // 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 (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ());
+  
+  // 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.
+  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
+  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
+  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);
+  std::vector <Point <dim> > current_points (n_components, Point <dim > ());
+  for (unsigned int support_point = 0; support_point < n_support_points; support_point++)
+    {
+      // setup valid data in the empty vectors
+      unsigned int component = 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
+  for (; cell != endc; cell++)
+    {
+      fe_values.reinit (cell);
+
+      for (unsigned int support_point = 0; support_point < n_support_points; support_point++)
+        {
+          unsigned int component = dof_handler->get_fe ().system_to_component_index (support_point).first;
+          Point<dim> test_point = fe_values.quadrature_point (support_point);
+
+          if (location.distance (test_point) < location.distance (current_points [component]))
+            {
+              // save the data
+              current_points [component] = test_point;
+              current_cell = cell;
+              current_fe_index [component] = support_point;
+            }
+        }
+    }
+
+
+  std::vector<unsigned int> 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.
+
+
+  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]]);
+    }
+
+  internal::PointValueHistory::PointGeometryData<dim> new_point_geometry_data (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
+      for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+        {
+          data_store_begin->second.push_back (std::vector<double> (0));
+        }
+    }
+}
+
+
+
+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!
+
+  // 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 (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ());
+  
+  // 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.
+  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
+  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
+
+  // 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);
+  std::vector <Point <dim> > temp_points (n_components, Point <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
+      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;
+    }
+  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
+  for (; cell != endc; cell++)
+    {
+      fe_values.reinit (cell);
+      for (unsigned int support_point = 0; support_point < n_support_points; support_point++)
+        {
+          unsigned int component = dof_handler->get_fe ().system_to_component_index (support_point).first;
+          Point<dim> test_point = fe_values.quadrature_point (support_point);
+
+          for (unsigned int point = 0; point < locations.size (); point++)
+            {
+              if (locations[point].distance (test_point) < locations[point].distance (current_points[point][component]))
+                {
+                  // save the data
+                  current_points[point][component] = test_point;
+                  current_cell[point] = cell;
+                  current_fe_index[point][component] = support_point;
+                }
+            }
+        }
+    }
+
+  std::vector<unsigned int> local_dof_indices (dof_handler->get_fe ().dofs_per_cell);
+  for (unsigned int point = 0; point < locations.size (); point++)
+    {
+      current_cell[point]->get_dof_indices (local_dof_indices);
+      std::vector <int> new_solution_indices;
+
+      for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+        {
+          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);
+
+      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
+          for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+            {
+              data_store_begin->second.push_back (std::vector<double> (0));
+            }
+        }
+    }
+}
+
+
+
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::add_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 () * (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);
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::close ()
+{
+  closed = true;
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::clear ()
+{
+  cleared = true;
+  dof_handler = SmartPointer<const DoFHandler<dim> > (0, "Cleared");
+}
+
+// Need to test that the internal data has a full and complete dataset for each key. That is that
+// the data has not got 'out of sync'. Testing that dataset_key is within 1 of independent_values
+// is cheap and is done in all three methods. Evaluate_field will check that its vector_name is within
+// 1 of dataset_key. However this leaves the possibility that the user has neglected to call evaluate_field
+// on one vector_name consistently. To catch this case start_new_dataset will call bool deap_check () which 
+// will test all vector_names and return a bool. This can be called from an Assert statement.
+
+
+
+template <int dim>
+template <class VECTOR>
+void PointValueHistory<dim>
+::evaluate_field (const std::string &vector_name, const VECTOR & solution)
+{
+  // 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 ()));
+  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 ());
+    }
+  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++)
+    {
+      // 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));
+        }
+    }
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::start_new_dataset (double key)
+{
+  // must be closed to add data to internal members.
+  Assert (closed, ExcInvalidState ());
+  Assert (!cleared, ExcInvalidState ());
+  Assert (deep_check (false), ExcDataLostSync ());
+
+  dataset_key.push_back (key);
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::push_back_independent (const std::vector <double> &indep_values)
+{
+  // 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));
+  Assert (n_indep != 0, ExcNoIndependent ());
+  Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
+
+  for (unsigned int component = 0; component < n_indep; component++)
+    independent_values[component].push_back (indep_values[component]);
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::write_gnuplot (const std::string &base_name)
+{
+  AssertThrow (closed, ExcInvalidState ());
+  AssertThrow (!cleared, ExcInvalidState ());
+  AssertThrow (deep_check (true), ExcDataLostSync ());
+
+  // write inputs to a file
+  if (n_indep != 0)
+    {
+      std::string filename = base_name + "_indep.gpl";
+      std::ofstream to_gnuplot (filename.c_str ());
+
+      to_gnuplot << "# Data independent of mesh location\n";
+
+      // write column headings
+      std::map <std::string, std::vector <std::vector <double> > >::iterator
+      data_store_begin = data_store.begin ();
+      to_gnuplot << "# <Key> ";
+
+      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++)
+        {
+          std::map <std::string, std::vector <std::vector <double> > >::iterator
+          data_store_begin = data_store.begin ();
+          to_gnuplot << dataset_key[key];
+
+          for (unsigned int component = 0; component < n_indep; component++)
+            {
+              to_gnuplot << " " << independent_values[component][key];
+            }
+          to_gnuplot << "\n";
+        }
+
+      to_gnuplot.close ();
+    }
+
+
+
+  // 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 ()));
+
+  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++)
+        {
+          // 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 upto 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";
+          for (unsigned int component = 0; component < n_components; component++)
+            {
+              to_gnuplot << "# " << node->solution_indices[component] << " : " << node->locations [component] << "\n";
+            }
+          to_gnuplot << "\n";
+
+          // write column headings
+          std::map <std::string, std::vector <std::vector <double> > >::iterator
+          data_store_begin = data_store.begin ();
+          to_gnuplot << "# <Key> ";
+
+          for (unsigned int component = 0; component < n_indep; component++)
+            {
+              to_gnuplot << "<Indep_" << component << "> ";
+            }
+
+          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 << "> ";
+                }
+            }
+          to_gnuplot << "\n";
+
+          // write data stored for the node
+          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 ();
+              to_gnuplot << dataset_key[key];
+
+              for (unsigned int component = 0; component < n_indep; component++)
+                {
+                  to_gnuplot << " " << independent_values[component][key];
+                }
+
+              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];
+                    }
+                }
+              to_gnuplot << "\n";
+            }
+
+          to_gnuplot.close ();
+        }
+    }
+}
+
+
+
+template <int dim>
+Vector<double> PointValueHistory<dim>
+::mark_locations ()
+{
+  // 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++)
+    {
+      for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+        {
+          dof_vector (node->solution_indices[component]) = 1;
+        }
+    }
+  return dof_vector;
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::get_points (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 ();
+
+  for (; node != point_geometry_data.end (); node++)
+    {
+      actual_points.push_back (node->locations);
+    }
+  locations = actual_points;
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::status (std::ostream &out)
+{
+  out << "***PointValueHistory status output***\n\n";
+  out << "Closed: " << closed << "\n";
+  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 ())
+    {
+      out << "No nodes 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";
+        }
+      }
+      else
+      {
+          out << "#Cannot access DoF_indices once cleared\n";
+      }
+    }
+  out << "\n";
+
+  std::cout << dataset_key.size () << " data keys are stored \n";
+
+
+  if (independent_values.size () != 0)
+    {
+      out << "Independent value(s): " << independent_values.size () << " : " << independent_values[0].size () << "\n";
+    }
+  else
+    {
+      out << "No independent values stored\n";
+    }
+
+  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++)
+    {
+      if (data_store_begin->second.size () != 0)
+        {
+          out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << (data_store_begin->second)[0].size () << "\n";
+        }
+      else
+        {
+          out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << "No points added" << "\n";
+        }
+    }
+  out << "\n";
+  out << "***end of status output***\n\n";
+}
+
+
+
+template <int dim>
+bool PointValueHistory<dim>
+::deep_check (bool strict)
+{
+  // test ways that it can fail, if control reaches last statement return true
+  if (strict)
+    {
+      if (n_indep != 0)
+        {
+          if (dataset_key.size () != independent_values[0].size ())
+            {
+              return false;
+            }
+        }
+      std::map <std::string, std::vector <std::vector <double> > >::iterator
+      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.
+            }
+        }
+      return true;
+    }
+  if (n_indep != 0)
+    {
+      if (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) >= 2)
+        {
+          return false;
+        }
+    }
+
+  if (n_dofs != 0)
+    {
+      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++)
+        {
+          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.
+        }
+    }
+  return true;
+}
+
+
+// explicit instantiations
+template class PointValueHistory<deal_II_dimension>;
+
+
+DEAL_II_NAMESPACE_CLOSE
+
diff --git a/deal.II/deal.II/source/numerics/point_value_history.inst.in b/deal.II/deal.II/source/numerics/point_value_history.inst.in
new file mode 100644 (file)
index 0000000..77fd01c
--- /dev/null
@@ -0,0 +1,7 @@
+for (VEC : SERIAL_VECTORS)
+{  
+  template
+  void PointValueHistory<deal_II_dimension>::evaluate_field 
+  (const std::string &,
+   const VEC &);
+}

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.