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