--- /dev/null
+//---------------------------------------------------------------------------
+//
+// 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 */
--- /dev/null
+//---------------------------------------------------------------------------
+//
+// 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
+