#include <deal.II/fe/fe_values.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/utilities.h>
+#include <numerics/data_postprocessor.h>
+#include <grid/grid_tools.h>
#include <vector>
#include <iostream>
class PointGeometryData
{
public:
- PointGeometryData(const std::vector <Point <dim> > &new_locations,
+ PointGeometryData(const
Point <dim> &new_requested_location, const std::vector <Point <dim> > &new_locations,
const std::vector <int> &new_sol_indices);
- std::vector <Point <dim> > locations;
+ Point <dim> requested_location;
+ std::vector <Point <dim> > support_point_locations;
std::vector <int> solution_indices;
};
}
* iterative process) graphs of solution values at specific points on the
* mesh. The user specifies the points which the solution should be monitored
* at ahead of time, as well as giving each solution vector that they want to
- * record a mnemonic name. Then for each step the user 'pushes back' the data
- * available from that time step and the class extracts data for the requested
- * points to store it. Finally once the computation is finished, the user can
+ * record a mnemonic name. Then, for each step the user calls one of the three
+ * available "evaluate field" methods to store the data from each time step,
+ * and the class extracts data for the requested
+ * points to store it. Finally, once the computation is finished, the user can
* request Gnuplot output files to be generated.
*
* The user can store extra variables which do not relate to mesh location by
* taken to solve the step and solver steps before convergence, or saving
* norms calculated.
*
- * Currently the code selects the nearest support points to a given point to
- * extract data from. This makes the code run at each time step very short,
+ * The three "evaluate field" methods each have different strengths and
+ * weaknesses making each suitable for different contexts:
+ * <ol>
+ * <li>Firstly, the @p evaluate_field version that does not take a @p DataPostprocessor object
+ * selects the nearest support points to a given point to
+ * extract data from. This makes the code that needs to be run at each time step very short,
* since looping over the mesh to extract the needed dof_index can be done
- * just once at the start. However this does lead to problems when
- * FiniteElements which do not assign dofs to actual mesh locations are used
- * (i.e. FEs without support points). The Functions::FEFieldFunction class
- * allows users to evaluate a solution at any point in a domain and will work
- * even for FEs without support points. This functionality is not currently
- * offered through this class since it has more overhead per iteration than
- * the approach implemented, but this class can be extended to offer it.
+ * just once at the start. However, this method is not suitable for
+ * FiniteElement objects that do not assign dofs to actual mesh locations
+ * (i.e. FEs without @ref GlossSupport "support points") or if adaptive mesh refinement is used.
+ * The class will throw an exception if the finite element does not have support points
+ * or any change in the @p dof_handler is identified. (The test is simply
+ * whether the number of dofs changes.)
+ *
+ * <li> Secondly, @p evaluate_field_at_requested_location calls @p
+ * VectorTools::point_value to compute values at the specific point requested.
+ * This method is valid for any FE that is supported by @p VectorTools::point_value.
+ * Specifically, this method can be called by codes using adaptive mesh refinement.
+ *
+ * <li>Finally, the class offers a function @p evaluate_field that takes a
+ * @p DataPostprocessor object. This method allows the deal.II
+ * data postprocessor to be used to compute new quantities from the
+ * solution on the fly. The values are located at the nearest quadrature
+ * point to the requested point. If the mesh is refined between calls, this
+ * point will change, so care must be taken when using this method in
+ * code using adaptive refinement, but as the output will be meaningful (in
+ * the sense that the quadrature point selected is guaranteed to remain in the
+ * same vicinity, the class does not prevent the use of this method in adaptive
+ * codes.
+ * </ol>
+ *
+ * When recording a new mnemonic name, the user must chose whether it should be
+ * a strictly scalar field (@p add_scalar_field_name) or whether it should
+ * depend on the @p FE (@p add_field_name). Calling @p add_field_name with a
+ * @p dof_handler using a vector @p FE will result in a vector field. The vector
+ * field will be of dimention @p FE.n_components(). Therefore, in this class the
+ * distinction between vector and scalar is as related to the @p FE and not
+ * whether a quantity is a logical vector.
*
* Following is a little code snippet that shows a common usage of this class:
*
* point is found and its details stored
* for use when @p evaluate_field is
* called. If more than one point is
- * required rather call this method as it
- * is more efficient than the add_points
+ * required, rather call this method as it
+ * is more efficient than the add_point
* method since it minimizes iterations
* over the mesh. The points are added to
* the internal database in the order
/**
* Put another mnemonic string (and hence
- * @p VECTOR) into the class. This also
- * adds extra entries for points that are
- * already in the class, so @p
- * add_field_name and @p add_points can
+ * @p VECTOR) into the class. This method
+ * adds storage space for FE.n_components()
+ * variables and should be used with methods
+ * that extract all FE components from a vector.
+ * This also adds extra entries for points
+ * that are already in the class, so
+ * @p add_field_name and @p add_points can
* be called in any order.
*/
void add_field_name(const std::string &vector_name);
+ /**
+ * Put another mnemonic string (and hence
+ * @p VECTOR) into the class. This method
+ * only adds storage space for a single
+ * variable whether the FE is scalar or not
+ * so it can be used with
+ * @p DataPostprocessor objects that compute
+ * a scalar value(s) from a vector field.
+ * This also adds extra entries for points
+ * that are already in the class, so
+ * @p add_field_name and @p add_points can
+ * be called in any order.
+ */
+ void add_scalar_field_name(const std::string &vector_name);
+
/**
* Extract values at the stored points
* from the VECTOR supplied and add them
- * to the new dataset in vector_name. If
- * a @p DoFHandler is used, this method
+ * to the new dataset in vector_name.
+ * The scalar_component input is only
+ * used if the FE has multiple components
+ * and the field to be evaluated is scalar.
+ * In this case, the scalar_component input
+ * is used to select the component used. If
+ * a @p DoFHandler is used, one (and only
+ * one) evaluate_field method
+ * must be called for each dataset (time
+ * step, iteration, etc) for each
+ * vector_name, otherwise a @p
+ * ExcDataLostSync error can occur.
+ */
+ template <class VECTOR>
+ void evaluate_field(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0);
+
+
+ /**
+ * Compute values using a DataPostprocessor object
+ * with the VECTOR supplied and add them
+ * to the new dataset in vector_name.
+ * The DataPostprocessor component_interpretation
+ * is used to determine select scalar or vector data.
+ * If the component_is_part_of_vector flag is used, then
+ * the indicated vector must have the same number of
+ * components as the FE associated with the dof_handler.
+ * The DataPostprocessor object is queried for field names.
+ * In the case of vector fields, the name of the first
+ * component is used as the overall name of the vector.
+ * The quadrature object supplied is used for all
+ * components of a vector field.
+ * If a @p DoFHandler is used, one (and only
+ * one) evaluate_field method
+ * must be called for each dataset (time
+ * step, iteration, etc) for each
+ * vector_name, otherwise a @p
+ * ExcDataLostSync error can occur.
+ */
+ template <class VECTOR>
+ void evaluate_field(const VECTOR & solution, const DataPostprocessor<dim> & data_postprocessor, const Quadrature<dim> & quadrature);
+
+
+ /**
+ * Extract values at the points actually
+ * requested from the VECTOR supplied and
+ * add them to the new dataset in vector_name.
+ * Unlike the other evaluate_field methods
+ * this method does not care if the dof_handler
+ * has been modified because it uses calls to
+ * @p VectorTools::point_value to extract there
+ * data. Therefore, if only this method is used,
+ * the class can be used with adaptive refinement.
+ * Specifically, the vector form of
+ * @p VectorTools::point_value that assumes a
+ * Q1 mapping is used.
+ * The scalar_component input is only
+ * used if the FE has multiple components
+ * and the field to be evaluated is scalar.
+ * In this case, the scalar_component input
+ * is used to select the component used. If
+ * a @p DoFHandler is used, one (and only
+ * one) evaluate_field method
* must be called for each dataset (time
* step, iteration, etc) for each
* vector_name, otherwise a @p
* ExcDataLostSync error can occur.
*/
template <class VECTOR>
- void evaluate_field(const std::string &vector_name, const VECTOR & solution);
+ void evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0);
/**
* Add the key for the current dataset to
* key and independent data. The file
* name agrees with the order the points
* were added to the class.
+ * The support point information is only
+ * meaningful if the dof_handler has not
+ * been changed. Therefore, if adaptive
+ * mesh refinement has been used the
+ * support point data should not be used.
+ * The optional parameter
+ * postprocessor_locations is used to
+ * add the postprocessor locations to the
+ * output files. If this is desired, the
+ * data should be obtained from a call to
+ * get_postprocessor_locations while the
+ * dof_handler is usable. The default
+ * parameter is an empty vector of strings,
+ * and will suppress postprocessor
+ * locations output.
*/
- void write_gnuplot(const std::string &base_name);
+ void write_gnuplot(const std::string &base_name, const std::vector <Point <dim> > postprocessor_locations = std::vector <Point <dim> > ());
/**
* data_out.attach_dof_handler(dof_handler);
*
* // Call the mark_locations method to get the vector with indices flagged
- * Vector<double> node_locations = node_monitor.mark_locations();
+ * Vector<double> support_point_locations = node_monitor.mark_locations();
*
* // Add the vector to the data_out object and write out a file in the usual way
- * data_out.add_data_vector(node_locations, "Monitor_Locations");
+ * data_out.add_data_vector(support_point_locations, "Monitor_Locations");
* data_out.build_patches(2);
* std::ofstream output("locations.gpl");
* data_out.write_gnuplot(output);
* @endcode
*/
+ Vector<double> mark_support_locations();
+
+
+ /**
+ * Depreciated:
+ *
+ * This function only exists for backward
+ * compatibility as this is the interface
+ * provided by previous versions of the library.
+ * The function mark_support_locations replaces
+ * it and reflects the fact that the locations
+ * marked are actually the support points.
+ */
Vector<double> mark_locations();
* the correct number of points by the
* method.
*/
+ void get_support_locations (std::vector <std::vector<Point <dim> > > & locations);
+
+ /**
+ * Depreciated:
+ *
+ * This function only exists for backward
+ * compatibility as this is the interface
+ * provided by previous versions of the library.
+ * The function get_support_locations replaces
+ * it and reflects the fact that the points
+ * returned are actually the support points.
+ */
void get_points (std::vector <std::vector<Point <dim> > > & locations);
+ /**
+ * Stores the actual location of the
+ * points used by the data_postprocessor.
+ * This can be used to compare with the
+ * points requested, for example by using
+ * the @p Point<dim>::distance function.
+ * Unlike the support_locations, these
+ * locations are computed every time the
+ * evaluate_field method is called with a
+ * postprocessor. This method uses the
+ * same algorithm so can will find the
+ * same points.
+ * For convenience, location is resized
+ * to the correct number of points by the
+ * method.
+ */
+ void get_postprocessor_locations (std::vector<Point <dim> > & locations, const Quadrature<dim> & quadrature);
/**
* Once datasets have been added to the
/**
* Saves data for each mnemonic entry.
- * data_store: mnemonic -> [component]
- * [key]
+ * data_store: mnemonic ->
+ * [point_0_components point_1_components
+ * ... point_n-1_components][key]
+ * This format facilitates scalar mnemonics
+ * in a vector space, because scalar mnemonics
+ * will only have one component per point.
+ * Vector components are strictly
+ * FE.n_components () long.
*/
std::map <std::string, std::vector <std::vector <double> > > data_store;
+ /**
+ * Saves a boolean representing whether
+ * the mnemonic links to a scalar field
+ * (true) or not (false).
+ */
+ std::map <std::string, bool> scalar_field;
+
/**
* Saves the location and other mesh info
* about support points.
*/
std::pair<std::string, std::vector <std::vector <double> > > pair_data; // could possibly be removed
- /**
- * Used to enforce @p closed state for some
- * methods.
- */
+ /**
+ * Used to enforce @p closed state for some
+ * methods.
+ */
bool closed;
- /**
- * Used to enforce @p !cleared state for
- * some methods.
- */
+ /**
+ * Used to enforce @p !cleared state for
+ * some methods.
+ */
bool cleared;
- /**
- * A smart pointer to the dof_handler
- * supplied to the constructor. This can be
- * released by calling @p clear().
- */
+ /**
+ * A smart pointer to the dof_handler
+ * supplied to the constructor. This can be
+ * released by calling @p clear().
+ */
SmartPointer<const DoFHandler<dim>,PointValueHistory<dim> > dof_handler;
/**
* Variable to check that number of dofs
* in the dof_handler does not change.
- */
- /* A cheap way to check that the
+ *
+ * A cheap way to check that the
* triangulation has not been refined in
* anyway since the class was set up.
* Refinement will invalidate stored dof
* indices linked to support points.
*/
unsigned int n_dofs;
- /**
- * Stores the number of independent
- * variables requested.
- */
+ /**
+ * Stores the number of independent
+ * variables requested.
+ */
unsigned int n_indep;
-
};
// $Id$
// Version: $Name$
//
-// Copyright (C) 2009, 2010, 2012 by Michael Rapson and the deal.II authors
+// Copyright (C) 2009, 2010 by Michael Rapson and the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_block_vector.h>
+#include <deal.II/numerics/vectors.h>
+
#include <deal.II/numerics/point_value_history.h>
/// Only a constructor needed for this class (a struct really)
template <int dim>
PointGeometryData<dim>
- ::PointGeometryData (const std::vector <Point <dim> > &new_locations,
- const std::vector <int> &new_sol_indices)
+ ::PointGeometryData (const Point <dim> &new_requested_location,
+ const std::vector <Point <dim> > &new_locations,
+ const std::vector <int> &new_sol_indices)
{
- locations = new_locations;
+ requested_location = new_requested_location;
+ support_point_locations = new_locations;
solution_indices = new_sol_indices;
}
}
template <int dim>
PointValueHistory<dim>
::PointValueHistory (const unsigned int n_independent_variables) :
- n_dofs (0),
- n_indep (n_independent_variables)
+ n_dofs (0),
+ n_indep (n_independent_variables)
{
closed = false;
cleared = false;
- // make a vector for keys
+ // make a vector for keys
dataset_key = std::vector <double> (); // initialize the std::vector
- // make a vector of independent values
+ // make a vector of independent values
independent_values
= std::vector<std::vector <double> > (n_indep, std::vector <double> (0));
}
template <int dim>
PointValueHistory<dim>::PointValueHistory (const DoFHandler<dim> & dof_handler,
- const unsigned int n_independent_variables) :
- dof_handler (&dof_handler),
- n_dofs (dof_handler.n_dofs ()),
- n_indep (n_independent_variables)
+ const unsigned int n_independent_variables) :
+ dof_handler (&dof_handler),
+ n_dofs (dof_handler.n_dofs ()),
+ n_indep (n_independent_variables)
{
closed = false;
cleared = false;
- // make a vector to store keys
+ // make a vector to store keys
dataset_key = std::vector <double> (); // initialize the std::vector
- // make a vector for the independent values
+ // make a vector for the independent values
independent_values
= std::vector<std::vector <double> > (n_indep, std::vector <double> (0));
}
dataset_key = point_value_history.dataset_key;
independent_values = point_value_history.independent_values;
data_store = point_value_history.data_store;
+ scalar_field = point_value_history.scalar_field;
point_geometry_data = point_value_history.point_geometry_data;
pair_data = point_value_history.pair_data;
closed = point_value_history.closed;
dataset_key = point_value_history.dataset_key;
independent_values = point_value_history.independent_values;
data_store = point_value_history.data_store;
+ scalar_field = point_value_history.scalar_field;
point_geometry_data = point_value_history.point_geometry_data;
pair_data = point_value_history.pair_data;
closed = point_value_history.closed;
void PointValueHistory<dim>
::add_point (const Point <dim> & location)
{
- // can't be closed to add additional points
- // or vectors
+ // can't be closed to add additional points
+ // or vectors
AssertThrow (!closed, ExcInvalidState ());
AssertThrow (!cleared, ExcInvalidState ());
AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
AssertThrow (n_dofs == dof_handler->n_dofs (),
- ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+ ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
- // Implementation assumes that support
- // points locations are dofs locations
+ // Implementation assumes that support
+ // points locations are dofs locations
AssertThrow (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ());
- // FEValues object to extract quadrature
- // points from
+ // FEValues object to extract quadrature
+ // points from
std::vector <Point <dim> >
unit_support_points = dof_handler->get_fe ().get_unit_support_points ();
- // While in general quadrature points seems
- // to refer to Gauss quadrature points, in
- // this case the quadrature points are
- // forced to be the support points of the
- // FE.
+ // While in general quadrature points seems
+ // to refer to Gauss quadrature points, in
+ // this case the quadrature points are
+ // forced to be the support points of the
+ // FE.
Quadrature<dim>
support_point_quadrature (dof_handler->get_fe ().get_unit_support_points ());
FEValues<dim> fe_values (dof_handler->get_fe (),
- support_point_quadrature,
- update_quadrature_points);
+ support_point_quadrature,
+ update_quadrature_points);
unsigned int n_support_points
= dof_handler->get_fe ().get_unit_support_points ().size ();
unsigned int n_components
= dof_handler->get_fe ().n_components ();
- // set up a loop over all the cells in the
- // DoFHandler
+ // set up a loop over all the cells in the
+ // DoFHandler
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler->begin_active ();
typename DoFHandler<dim>::active_cell_iterator
endc = dof_handler->end ();
- // default values to be replaced as closer
- // points are found however they need to be
- // consistent in case they are actually
- // chosen
+ // default values to be replaced as closer
+ // points are found however they need to be
+ // consistent in case they are actually
+ // chosen
typename DoFHandler<dim>::active_cell_iterator current_cell = cell;
std::vector <unsigned int> current_fe_index (n_components, 0); // need one index per component
fe_values.reinit (cell);
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>
- new_point_geometry_data (current_points, new_solution_indices);
+ new_point_geometry_data (location, current_points, new_solution_indices);
point_geometry_data.push_back (new_point_geometry_data);
std::map <std::string, std::vector <std::vector <double> > >::iterator
data_store_begin = data_store.begin ();
for (; data_store_begin != data_store.end (); data_store_begin++)
{
- // add an extra row to each vector
- // entry
+ // add an extra row to each vector
+ // entry
+ bool scalar = (scalar_field.find (data_store_begin->first))->second;
for (unsigned int component = 0;
- component < dof_handler->get_fe ().n_components (); component++)
+ component < (scalar ? 1 : dof_handler->get_fe ().n_components ()); component++)
{
data_store_begin->second.push_back (std::vector<double> (0));
}
void PointValueHistory<dim>
::add_points (const std::vector <Point <dim> > & locations)
{
- // This algorithm adds points in the same
- // order as they appear in the vector
- // locations and users may depend on this
- // so do not change order added!
+ // This algorithm adds points in the same
+ // order as they appear in the vector
+ // locations and users may depend on this
+ // so do not change order added!
- // can't be closed to add additional points or vectors
+ // can't be closed to add additional points or vectors
AssertThrow (!closed, ExcInvalidState ());
AssertThrow (!cleared, ExcInvalidState ());
AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
- // Implementation assumes that support
- // points locations are dofs locations
+ // Implementation assumes that support
+ // points locations are dofs locations
AssertThrow (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ());
- // FEValues object to extract quadrature
- // points from
+ // FEValues object to extract quadrature
+ // points from
std::vector <Point <dim> > unit_support_points = dof_handler->get_fe ().get_unit_support_points ();
- // While in general quadrature points seems
- // to refer to Gauss quadrature points, in
- // this case the quadrature points are
- // forced to be the support points of the
- // FE.
+ // While in general quadrature points seems
+ // to refer to Gauss quadrature points, in
+ // this case the quadrature points are
+ // forced to be the support points of the
+ // FE.
Quadrature<dim> support_point_quadrature (dof_handler->get_fe ().get_unit_support_points ());
FEValues<dim> fe_values (dof_handler->get_fe (), support_point_quadrature, update_quadrature_points);
unsigned int n_support_points = dof_handler->get_fe ().get_unit_support_points ().size ();
unsigned int n_components = dof_handler->get_fe ().n_components ();
- // set up a loop over all the cells in the
- // DoFHandler
+ // set up a loop over all the cells in the
+ // DoFHandler
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler->begin_active ();
typename DoFHandler<dim>::active_cell_iterator endc = dof_handler->end ();
- // default values to be replaced as closer
- // points are found however they need to be
- // consistent in case they are actually
- // chosen vector <vector>s defined where
- // previously single vectors were used
+ // default values to be replaced as closer
+ // points are found however they need to be
+ // consistent in case they are actually
+ // chosen vector <vector>s defined where
+ // previously single vectors were used
- // need to store one value per point per component
+ // need to store one value per point per component
std::vector <typename DoFHandler<dim>::active_cell_iterator > current_cell (locations.size (), cell);
fe_values.reinit (cell);
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;
new_solution_indices.push_back (local_dof_indices[current_fe_index[point][component]]);
}
- internal::PointValueHistory::PointGeometryData<dim> new_point_geometry_data (current_points[point], new_solution_indices);
+ internal::PointValueHistory::PointGeometryData<dim> new_point_geometry_data (locations[point], current_points[point], new_solution_indices);
point_geometry_data.push_back (new_point_geometry_data);
std::map <std::string, std::vector <std::vector <double> > >::iterator
- data_store_begin = data_store.begin ();
+ data_store_begin = data_store.begin ();
for (; data_store_begin != data_store.end (); data_store_begin++)
{
- // add an extra row to each vector
- // entry
- for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+ // add an extra row to each vector
+ // entry
+ bool scalar = (scalar_field.find (data_store_begin->first))->second;
+ for (unsigned int component = 0; component < (scalar ? 1 : dof_handler->get_fe ().n_components ()); component++)
{
data_store_begin->second.push_back (std::vector<double> (0));
}
void PointValueHistory<dim>
::add_field_name (const std::string &vector_name)
{
- // can't be closed to add additional points
- // or vectors
+ // can't be closed to add additional points
+ // or vectors
AssertThrow (!closed, ExcInvalidState ());
AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
AssertThrow (!cleared, ExcInvalidState ());
AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
- // make and add a new vector
- // point_geometry_data.size() long
+ // make and add a new vector
+ // point_geometry_data.size() long
pair_data.first = vector_name;
int n_datastreams = point_geometry_data.size () * (dof_handler->get_fe ().n_components ()); // each point has n_components sub parts
std::vector < std::vector <double> > vector_size (n_datastreams,
std::vector <double> (0));
pair_data.second = vector_size;
data_store.insert (pair_data);
+ // this is only a scalar field if
+ // the FE is scalar
+ if (dof_handler->get_fe ().n_components () == 1)
+ scalar_field.insert (std::pair <std::string, bool> (vector_name, true));
+ else
+ scalar_field.insert (std::pair <std::string, bool> (vector_name, false));
+}
+
+
+
+template <int dim>
+void PointValueHistory<dim>
+::add_scalar_field_name (const std::string &vector_name)
+{
+ // can't be closed to add additional points
+ // or vectors
+ AssertThrow (!closed, ExcInvalidState ());
+ AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
+ AssertThrow (!cleared, ExcInvalidState ());
+ AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+
+
+ // make and add a new vector
+ // point_geometry_data.size() long
+ pair_data.first = vector_name;
+ int n_datastreams = point_geometry_data.size (); // each point has 1 value for this function
+ std::vector < std::vector <double> > vector_size (n_datastreams,
+ std::vector <double> (0));
+ pair_data.second = vector_size;
+ data_store.insert (pair_data);
+ scalar_field.insert (std::pair <std::string, bool> (vector_name, true)); // scalar by definition in this method
}
template <int dim>
template <class VECTOR>
void PointValueHistory<dim>
-::evaluate_field (const std::string &vector_name, const VECTOR & solution)
+::evaluate_field (const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component)
{
- // must be closed to add data to internal
- // members.
+ // must be closed to add data to internal
+ // members.
Assert (closed, ExcInvalidState ());
Assert (!cleared, ExcInvalidState ());
Assert (n_dofs != 0, ExcDoFHandlerRequired ());
{
Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
}
- typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
- for (unsigned int data_store_index = 0; node != point_geometry_data.end (); node++, data_store_index++)
+ // Look up the field name and get an
+ // iterator for the map. Doing this
+ // up front means that it only needs
+ // to be done once and also allows us
+ // to check vector_name is in the map.
+ typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(vector_name);
+ Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+ for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
{
- // step through each node, to access
- // the Solution_index and the vector
- // index
- for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
- {
- unsigned int solution_index = node->solution_indices[component];
- (data_store[vector_name])[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (solution (solution_index));
- }
+ if (!(scalar_field.find (vector_name))->second)
+ {
+ // The solution_indices data has a
+ // one-to-one correspondance to FE
+ // components therefore we simply
+ // step through each point, to access
+ // the Solution_index and the vector
+ // index
+ for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+ {
+ unsigned int solution_index = point->solution_indices[component];
+ data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (solution (solution_index));
+ }
+ }
+ else
+ {
+ // We have found a scalar field
+ // added to a dof_handler with a
+ // vector FE. Therefore, we take
+ // the value specified by the user
+ // in the scalar_component input.
+ Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError());
+ unsigned int solution_index = point->solution_indices[scalar_component];
+ data_store_field->second[data_store_index].push_back (solution (solution_index));
+ }
}
}
+
+
+template <int dim>
+template <class VECTOR>
+void PointValueHistory<dim>
+::evaluate_field(const VECTOR & solution, const DataPostprocessor< dim> & data_postprocessor, const Quadrature<dim> & quadrature)
+{
+ // must be closed to add data to internal
+ // members.
+ Assert (closed, ExcInvalidState ());
+ Assert (!cleared, ExcInvalidState ());
+ Assert (n_dofs != 0, ExcDoFHandlerRequired ());
+ if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0
+ {
+ Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
+ }
+
+// Make an FEValues object
+ const UpdateFlags update_flags = data_postprocessor.get_needed_update_flags() | update_quadrature_points;
+ Assert (!(update_flags & update_normal_vectors),
+ ExcMessage("The update of normal vectors may not be requested for evaluation of "
+ "data on cells via DataPostprocessor."));
+ FEValues<dim> fe_values (dof_handler->get_fe (), quadrature, update_flags);
+ unsigned int n_components = dof_handler->get_fe ().n_components ();
+ unsigned int n_quadrature_points = quadrature.size();
+
+ // extract the data names, used to size the data_postprocessor output
+ std::vector<std::string> field_names = data_postprocessor.get_names ();
+
+// Loop over points and find correct cell
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+ for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
+ {
+ // we now have a point to query, need to know what cell it is in
+ Point <dim> requested_location = point->requested_location;
+ typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1<dim>(), *dof_handler, requested_location).first;
+
+
+ fe_values.reinit (cell);
+ std::vector< Vector< double > > computed_quantities (1, Vector <double> (field_names.size())); // just one point needed
+
+ // The case of a scalar FE
+ if (n_components == 1)
+ {
+ // Extract data for the PostProcessor object
+ std::vector< double > uh (n_quadrature_points, 0.0);
+ std::vector< Tensor< 1, dim > > duh (n_quadrature_points, Tensor <1, dim> ());
+ std::vector< Tensor< 2, dim > > dduh (n_quadrature_points, Tensor <2, dim> ());
+ std::vector<Point<dim> > dummy_normals (1, Point<dim> ());
+ std::vector<Point<dim> > evaluation_points;
+ // at each point there is
+ // only one component of
+ // value, gradient etc.
+ if (update_flags & update_values)
+ fe_values.get_function_values (solution,
+ uh);
+ if (update_flags & update_gradients)
+ fe_values.get_function_gradients (solution,
+ duh);
+ if (update_flags & update_hessians)
+ fe_values.get_function_hessians (solution,
+ dduh);
+
+ // find the closest quadrature point
+ evaluation_points = fe_values.get_quadrature_points();
+ double distance = cell->diameter ();
+ unsigned int selected_point = 0;
+ for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
+ {
+ if (requested_location.distance (evaluation_points[q_point]) < distance)
+ {
+ selected_point = q_point;
+ distance = requested_location.distance (evaluation_points[q_point]);
+ }
+ }
+
+ // Call compute_derived_quantities_vector or compute_derived_quantities_scalar
+ data_postprocessor.
+ compute_derived_quantities_scalar(std::vector< double > (1, uh[selected_point]),
+ std::vector< Tensor< 1, dim > > (1, duh[selected_point]),
+ std::vector< Tensor< 2, dim > > (1, dduh[selected_point]),
+ dummy_normals,
+ std::vector<Point<dim> > (1, evaluation_points[selected_point]),
+ computed_quantities);
+
+ }
+ else // The case of a vector FE
+ {
+ // Extract data for the PostProcessor object
+ std::vector< Vector< double > > uh (n_quadrature_points, Vector <double> (n_components));
+ std::vector< std::vector< Tensor< 1, dim > > > duh (n_quadrature_points, std::vector< Tensor< 1, dim > > (n_components, Tensor< 1, dim >()));
+ std::vector< std::vector< Tensor< 2, dim > > > dduh (n_quadrature_points, std::vector< Tensor< 2, dim > > (n_components, Tensor< 2, dim >()));
+ std::vector<Point<dim> > dummy_normals (1, Point<dim> ());
+ std::vector<Point<dim> > evaluation_points;
+ // at each point there is
+ // a vector valued
+ // function and its
+ // derivative...
+ if (update_flags & update_values)
+ fe_values.get_function_values (solution,
+ uh);
+ if (update_flags & update_gradients)
+ fe_values.get_function_gradients (solution,
+ duh);
+ if (update_flags & update_hessians)
+ fe_values.get_function_hessians (solution,
+ dduh);
+
+ // find the closest quadrature point
+ evaluation_points = fe_values.get_quadrature_points();
+ double distance = cell->diameter ();
+ unsigned int selected_point = 0;
+ for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
+ {
+ if (requested_location.distance (evaluation_points[q_point]) < distance)
+ {
+ selected_point = q_point;
+ distance = requested_location.distance (evaluation_points[q_point]);
+ }
+ }
+
+ // Call compute_derived_quantities_vector or compute_derived_quantities_scalar
+ data_postprocessor.
+ compute_derived_quantities_vector(std::vector< Vector< double > > (1, uh[selected_point]),
+ std::vector< std::vector< Tensor< 1, dim > > > (1, duh[selected_point]),
+ std::vector< std::vector< Tensor< 2, dim > > > (1, dduh[selected_point]),
+ dummy_normals,
+ std::vector<Point<dim> > (1, evaluation_points[selected_point]),
+ computed_quantities);
+ }
+
+ // we now have the data and need to save it
+// std::vector<DataComponentInterpretation::DataComponentInterpretation> component_interpretation = data_postprocessor.get_data_component_interpretation ();
+
+ // loop over data names
+ typename std::vector<std::string>::iterator name = field_names.begin();
+ for (unsigned int computed_quantities_index = 0; name != field_names.end();)
+ {
+ // this loop is incremented in
+ // separate locations depending on
+ // whether the name is a vector_field
+ // scalar_field or not in the
+ // data_store at all.
+ typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(*name);
+ if (data_store_field != data_store.end()) // Ignore names not in the class
+ {
+ if (!(scalar_field.find (*name))->second)
+ {
+ // The solution_indices data has a
+ // one-to-one correspondance to FE
+ // components therefore we simply
+ // step through each point, to access
+ // the Solution_index and the vector
+ // index
+ for (unsigned int component = 0; component < n_components; component++, name++, computed_quantities_index++)
+ {
+ Assert (computed_quantities_index < field_names.size(), ExcMessage ("Insufficient components in DataPostprocessor for current vector field."));
+
+ data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (computed_quantities[0](computed_quantities_index));
+ }
+ }
+ else
+ {
+ // We have found a scalar field
+ // added to a dof_handler with a
+ // vector FE. Therefore, we take
+ // the value specified by the user
+ // in the scalar_component input.
+ Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError());
+ data_store_field->second[data_store_index].push_back (computed_quantities[0](computed_quantities_index));
+ name ++;
+ computed_quantities_index ++;
+ }
+ }
+ else
+ {
+ name ++;
+ computed_quantities_index ++;
+ }
+ }
+ } // end of loop over points
+}
+
+
+
+
+
+
+template <int dim>
+template <class VECTOR>
+void PointValueHistory<dim>
+::evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component)
+{
+ // must be closed to add data to internal
+ // members.
+ Assert (closed, ExcInvalidState ());
+ Assert (!cleared, ExcInvalidState ());
+ Assert (n_dofs != 0, ExcDoFHandlerRequired ());
+ if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0
+ {
+ Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
+ }
+ // Look up the field name and get an
+ // iterator for the map. Doing this
+ // up front means that it only needs
+ // to be done once and also allows us
+ // to check vector_name is in the map.
+ typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(vector_name);
+ Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+ Vector <double> value (dof_handler->get_fe().n_components());
+ for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
+ {
+ // Make a Vector <double> for the value
+ // at the point. It will have as many
+ // components as there are in the fe.
+ VectorTools::point_value (*dof_handler, solution, point->requested_location, value);
+
+ if (!(scalar_field.find (vector_name))->second)
+ {
+ // The solution_indices data has a
+ // one-to-one correspondance to FE
+ // components therefore we simply
+ // step through each point, to access
+ // the Solution_index and the vector
+ // index
+ for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+ {
+ data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (value (component));
+ }
+ }
+ else
+ {
+ // We have found a scalar field
+ // added to a dof_handler with a
+ // vector FE. Therefore, we take
+ // the value specified by the user
+ // in the scalar_component input.
+ Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError());
+ data_store_field->second[data_store_index].push_back (value (scalar_component));
+ }
+ }
+}
+
+
template <int dim>
void PointValueHistory<dim>
::start_new_dataset (double key)
{
- // must be closed to add data to internal
- // members.
+ // must be closed to add data to internal
+ // members.
Assert (closed, ExcInvalidState ());
Assert (!cleared, ExcInvalidState ());
Assert (deep_check (false), ExcDataLostSync ());
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));
template <int dim>
void PointValueHistory<dim>
-::write_gnuplot (const std::string &base_name)
+::write_gnuplot (const std::string &base_name, const std::vector <Point <dim> > postprocessor_locations)
{
AssertThrow (closed, ExcInvalidState ());
AssertThrow (!cleared, ExcInvalidState ());
AssertThrow (deep_check (true), ExcDataLostSync ());
- // write inputs to a file
+ // write inputs to a file
if (n_indep != 0)
{
std::string filename = base_name + "_indep.gpl";
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 << "\n";
- // write general data stored
+ // write general data stored
for (unsigned int key = 0; key < dataset_key.size (); key++)
{
to_gnuplot << dataset_key[key];
- // write points to a file
+ // write points to a file
if (n_dofs != 0)
{
AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
- AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
+ AssertThrow (postprocessor_locations.size() == 0 || postprocessor_locations.size() == point_geometry_data.size(), ExcDimensionMismatch (postprocessor_locations.size(), point_geometry_data.size()));
+ // We previously required the
+ // number of dofs to remain the
+ // same to provide some sort of
+ // test on the relevance of the
+ // support point indices stored.
+ // We now relax that to allow
+ // adaptive refinement strategies
+ // to make use of the
+ // evaluate_field_requested_locations
+ // method. Note that the support point
+ // information is not meaningful if
+ // the number of dofs has changed.
+ //AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
unsigned int n_components = dof_handler->get_fe ().n_components ();
- typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
- for (unsigned int data_store_index = 0; node != point_geometry_data.end (); node++, data_store_index++)
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+ for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
{
- // for each point, open a file to
- // be written to
+ // for each point, open a file to
+ // be written to
std::string filename = base_name + "_" + Utilities::int_to_string (data_store_index, 2) + ".gpl"; // store by order pushed back
- // due to
- // Utilities::int_to_string(data_store_index,
- // 2) call, can handle up to 100
- // points
+ // due to
+ // Utilities::int_to_string(data_store_index,
+ // 2) call, can handle up to 100
+ // points
std::ofstream to_gnuplot (filename.c_str ());
- // put helpful info about the
- // support point into the file as
- // comments
- to_gnuplot << "# DoF_index : Location (for each component)\n";
+ // put helpful info about the
+ // support point into the file as
+ // comments
+ to_gnuplot << "# Requested location: " << point->requested_location << "\n";
+ to_gnuplot << "# DoF_index : Support location (for each component)\n";
for (unsigned int component = 0; component < n_components; component++)
{
- to_gnuplot << "# " << node->solution_indices[component] << " : " << node->locations [component] << "\n";
+ to_gnuplot << "# " << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n";
}
- to_gnuplot << "\n";
+ if (postprocessor_locations.size() != 0)
+ {
+ to_gnuplot << "# Postprocessor location: " << postprocessor_locations[data_store_index] << "\n";
+ }
+ to_gnuplot << "#\n";
- // write column headings
+ // write column headings
std::map <std::string, std::vector <std::vector <double> > >::iterator
- data_store_begin = data_store.begin ();
+ data_store_begin = data_store.begin ();
to_gnuplot << "# <Key> ";
for (unsigned int component = 0; component < n_indep; component++)
for (; data_store_begin != data_store.end (); data_store_begin++)
{
- for (unsigned int component = 0; component < n_components; component++)
- {
- to_gnuplot << "<" << data_store_begin->first << "_" << component << "> ";
- }
+ if (!(scalar_field.find (data_store_begin->first))->second)
+ for (unsigned int component = 0; component < n_components; component++)
+ {
+ to_gnuplot << "<" << data_store_begin->first << "_" << component << "> ";
+ }
+ else
+ to_gnuplot << "<" << data_store_begin->first << "_" << 0 << "> ";
}
to_gnuplot << "\n";
- // write data stored for the node
+ // write data stored for the point
for (unsigned int key = 0; key < dataset_key.size (); key++)
{
std::map <std::string, std::vector <std::vector <double> > >::iterator
- data_store_begin = data_store.begin ();
+ data_store_begin = data_store.begin ();
to_gnuplot << dataset_key[key];
for (unsigned int component = 0; component < n_indep; component++)
for (; data_store_begin != data_store.end (); data_store_begin++)
{
- for (unsigned int component = 0; component < n_components; component++)
- {
- to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key];
- }
+ bool scalar = (scalar_field.find (data_store_begin->first))->second;
+ if (scalar)
+ {
+ to_gnuplot << " " << (data_store_begin->second)[data_store_index][key];
+ }
+ else
+ {
+ for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+ {
+ to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key];
+ }
+ }
}
to_gnuplot << "\n";
}
template <int dim>
Vector<double> PointValueHistory<dim>
-::mark_locations ()
+::mark_support_locations ()
{
- // a method to put a one at each point on
- // the grid where a location is defined
+ // a method to put a one at each point on
+ // the grid where a location is defined
AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
AssertThrow (!cleared, ExcInvalidState ());
AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
Vector<double> dof_vector (dof_handler->n_dofs ());
- typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
- for (; node != point_geometry_data.end (); node++)
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+ for (; point != point_geometry_data.end (); point++)
{
for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
{
- dof_vector (node->solution_indices[component]) = 1;
+ dof_vector (point->solution_indices[component]) = 1;
}
}
return dof_vector;
}
+template <int dim>
+Vector<double> PointValueHistory<dim>
+::mark_locations ()
+{
+ return mark_support_locations ();
+}
+
template <int dim>
void PointValueHistory<dim>
-::get_
points (std::vector <std::vector<Point <dim> > > & locations)
+::get_
support_locations (std::vector <std::vector<Point <dim> > > & locations)
{
AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ());
AssertThrow (!cleared, ExcInvalidState ());
AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ()));
std::vector <std::vector <Point <dim> > > actual_points;
- typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
- for (; node != point_geometry_data.end (); node++)
+ for (; point != point_geometry_data.end (); point++)
{
- actual_points.push_back (node->locations);
+ actual_points.push_back (point->support_point_locations);
}
locations = actual_points;
}
+template <int dim>
+void PointValueHistory<dim>
+::get_points (std::vector <std::vector<Point <dim> > > & locations)
+{
+ get_support_locations (locations);
+}
+
+
+template <int dim>
+void PointValueHistory<dim>
+::get_postprocessor_locations (std::vector<Point <dim> > & locations, const Quadrature<dim> & quadrature)
+{
+ Assert (!cleared, ExcInvalidState ());
+ Assert (n_dofs != 0, ExcDoFHandlerRequired ());
+
+ locations = std::vector<Point <dim> > ();
+
+ FEValues<dim> fe_values (dof_handler->get_fe (), quadrature, update_quadrature_points);
+ unsigned int n_quadrature_points = quadrature.size();
+ std::vector<Point<dim> > evaluation_points;
+
+// Loop over points and find correct cell
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+ for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
+ {
+ // we now have a point to query, need to know what cell it is in
+ Point <dim> requested_location = point->requested_location;
+ typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1<dim>(), *dof_handler, requested_location).first;
+ fe_values.reinit (cell);
+
+ evaluation_points = fe_values.get_quadrature_points();
+ double distance = cell->diameter ();
+ unsigned int selected_point = 0;
+
+ for (unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
+ {
+ if (requested_location.distance (evaluation_points[q_point]) < distance)
+ {
+ selected_point = q_point;
+ distance = requested_location.distance (evaluation_points[q_point]);
+ }
+ }
+
+ locations.push_back (evaluation_points[selected_point]);
+ }
+}
+
template <int dim>
void PointValueHistory<dim>
out << "Cleared: " << cleared << "\n";
out << "Geometric Data" << "\n";
- typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator node = point_geometry_data.begin ();
- if (node == point_geometry_data.end ())
+ typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
+ if (point == point_geometry_data.end ())
{
- out << "No nodes stored currently\n";
+ out << "No points stored currently\n";
}
else
{
if (!cleared)
- {
- out << "# DoF_index : Location (for each component)\n";
- for (; node != point_geometry_data.end (); node++)
- {
- for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
- {
- out << node->solution_indices[component] << " : " << node->locations [component] << "\n";
- }
- out << "\n";
- }
- }
+ {
+ for (; point != point_geometry_data.end (); point++)
+ {
+ out << "# Requested location: " << point->requested_location << "\n";
+ out << "# DoF_index : Support location (for each component)\n";
+ for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
+ {
+ out << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n";
+ }
+ out << "\n";
+ }
+ }
else
- {
+ {
out << "#Cannot access DoF_indices once cleared\n";
- }
+ }
}
out << "\n";
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 ();
+ data_store_begin = data_store.begin ();
if (n_dofs != 0)
{
for (; data_store_begin != data_store.end (); data_store_begin++)
{
if ((data_store_begin->second)[0].size () != dataset_key.size ())
return false;
- // this loop only tests one
- // member for each name,
- // i.e. checks the user it will
- // not catch internal errors
- // which do not update all
- // fields for a name.
+ // this loop only tests one
+ // member for each name,
+ // i.e. checks the user it will
+ // not catch internal errors
+ // which do not update all
+ // fields for a name.
}
}
return true;
if (n_dofs != 0)
{
std::map <std::string, std::vector <std::vector <double> > >::iterator
- data_store_begin = data_store.begin ();
+ data_store_begin = data_store.begin ();
for (; data_store_begin != data_store.end (); data_store_begin++)
{
if (std::abs ((int) (data_store_begin->second)[0].size () - (int) dataset_key.size ()) >= 2)
return false;
- // this loop only tests one member
- // for each name, i.e. checks the
- // user it will not catch internal
- // errors which do not update all
- // fields for a name.
+ // this loop only tests one member
+ // for each name, i.e. checks the
+ // user it will not catch internal
+ // errors which do not update all
+ // fields for a name.
}
}
return true;
projects / dealii-svn.git / commitdiff