From 286be4b26bc66343ba4b6ad84660218370500375 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 1 Jun 2012 19:07:07 +0000 Subject: [PATCH] Patch by Michael Rapson: Major update of the PointValueHistory class. git-svn-id: https://svn.dealii.org/trunk@25595 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/numerics/point_value_history.h | 275 ++++-- .../source/numerics/point_value_history.cc | 780 +++++++++++++----- .../numerics/point_value_history.inst.in | 28 +- 3 files changed, 837 insertions(+), 246 deletions(-) diff --git a/deal.II/include/deal.II/numerics/point_value_history.h b/deal.II/include/deal.II/numerics/point_value_history.h index 9ab521fe12..5d3b023cfa 100644 --- a/deal.II/include/deal.II/numerics/point_value_history.h +++ b/deal.II/include/deal.II/numerics/point_value_history.h @@ -26,6 +26,8 @@ #include #include #include +#include +#include #include #include @@ -59,9 +61,10 @@ namespace internal class PointGeometryData { public: - PointGeometryData(const std::vector > &new_locations, + PointGeometryData(const Point &new_requested_location, const std::vector > &new_locations, const std::vector &new_sol_indices); - std::vector > locations; + Point requested_location; + std::vector > support_point_locations; std::vector solution_indices; }; } @@ -74,9 +77,10 @@ namespace internal * iterative process) graphs of solution values at specific points on the * mesh. The user specifies the points which the solution should be monitored * at ahead of time, as well as giving each solution vector that they want to - * record a mnemonic name. Then for each step the user 'pushes back' the data - * available from that time step and the class extracts data for the requested - * points to store it. Finally once the computation is finished, the user can + * record a mnemonic name. Then, for each step the user calls one of the three + * available "evaluate field" methods to store the data from each time step, + * and the class extracts data for the requested + * points to store it. Finally, once the computation is finished, the user can * request Gnuplot output files to be generated. * * The user can store extra variables which do not relate to mesh location by @@ -87,16 +91,44 @@ namespace internal * taken to solve the step and solver steps before convergence, or saving * norms calculated. * - * Currently the code selects the nearest support points to a given point to - * extract data from. This makes the code run at each time step very short, + * The three "evaluate field" methods each have different strengths and + * weaknesses making each suitable for different contexts: + *
    + *
  1. 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.) + * + *
  2. 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. + * + *
  3. 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. + *
+ * + * When recording a new mnemonic name, the user must chose whether it should be + * a strictly scalar field (@p add_scalar_field_name) or whether it should + * depend on the @p FE (@p add_field_name). Calling @p add_field_name with a + * @p dof_handler using a vector @p FE will result in a vector field. The vector + * field will be of dimention @p FE.n_components(). Therefore, in this class the + * distinction between vector and scalar is as related to the @p FE and not + * whether a quantity is a logical vector. * * Following is a little code snippet that shows a common usage of this class: * @@ -236,8 +268,8 @@ class PointValueHistory * point is found and its details stored * for use when @p evaluate_field is * called. If more than one point is - * required rather call this method as it - * is more efficient than the add_points + * required, rather call this method as it + * is more efficient than the add_point * method since it minimizes iterations * over the mesh. The points are added to * the internal database in the order @@ -253,26 +285,104 @@ class PointValueHistory /** * Put another mnemonic string (and hence - * @p VECTOR) into the class. This also - * adds extra entries for points that are - * already in the class, so @p - * add_field_name and @p add_points can + * @p VECTOR) into the class. This method + * adds storage space for FE.n_components() + * variables and should be used with methods + * that extract all FE components from a vector. + * This also adds extra entries for points + * that are already in the class, so + * @p add_field_name and @p add_points can * be called in any order. */ void add_field_name(const std::string &vector_name); + /** + * Put another mnemonic string (and hence + * @p VECTOR) into the class. This method + * only adds storage space for a single + * variable whether the FE is scalar or not + * so it can be used with + * @p DataPostprocessor objects that compute + * a scalar value(s) from a vector field. + * This also adds extra entries for points + * that are already in the class, so + * @p add_field_name and @p add_points can + * be called in any order. + */ + void add_scalar_field_name(const std::string &vector_name); + /** * Extract values at the stored points * from the VECTOR supplied and add them - * to the new dataset in vector_name. If - * a @p DoFHandler is used, this method + * to the new dataset in vector_name. + * The scalar_component input is only + * used if the FE has multiple components + * and the field to be evaluated is scalar. + * In this case, the scalar_component input + * is used to select the component used. If + * a @p DoFHandler is used, one (and only + * one) evaluate_field method + * must be called for each dataset (time + * step, iteration, etc) for each + * vector_name, otherwise a @p + * ExcDataLostSync error can occur. + */ + template + 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 + void evaluate_field(const VECTOR & solution, const DataPostprocessor & data_postprocessor, const Quadrature & 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 - void evaluate_field(const std::string &vector_name, const VECTOR & solution); + void evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0); /** * Add the key for the current dataset to @@ -313,8 +423,23 @@ class PointValueHistory * key and independent data. The file * name agrees with the order the points * were added to the class. + * The support point information is only + * meaningful if the dof_handler has not + * been changed. Therefore, if adaptive + * mesh refinement has been used the + * support point data should not be used. + * The optional parameter + * postprocessor_locations is used to + * add the postprocessor locations to the + * output files. If this is desired, the + * data should be obtained from a call to + * get_postprocessor_locations while the + * dof_handler is usable. The default + * parameter is an empty vector of strings, + * and will suppress postprocessor + * locations output. */ - void write_gnuplot(const std::string &base_name); + void write_gnuplot(const std::string &base_name, const std::vector > postprocessor_locations = std::vector > ()); /** @@ -336,15 +461,28 @@ class PointValueHistory * data_out.attach_dof_handler(dof_handler); * * // Call the mark_locations method to get the vector with indices flagged - * Vector node_locations = node_monitor.mark_locations(); + * Vector 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 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 mark_locations(); @@ -359,8 +497,37 @@ class PointValueHistory * the correct number of points by the * method. */ + void get_support_locations (std::vector > > & 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 > > & 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::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 > & locations, const Quadrature & quadrature); /** * Once datasets have been added to the @@ -490,11 +657,24 @@ class PointValueHistory /** * Saves data for each mnemonic entry. - * data_store: mnemonic -> [component] - * [key] + * data_store: mnemonic -> + * [point_0_components point_1_components + * ... point_n-1_components][key] + * This format facilitates scalar mnemonics + * in a vector space, because scalar mnemonics + * will only have one component per point. + * Vector components are strictly + * FE.n_components () long. */ std::map > > data_store; + /** + * Saves a boolean representing whether + * the mnemonic links to a scalar field + * (true) or not (false). + */ + std::map scalar_field; + /** * Saves the location and other mesh info * about support points. @@ -508,43 +688,42 @@ class PointValueHistory */ std::pair > > 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,PointValueHistory > 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; - }; diff --git a/deal.II/source/numerics/point_value_history.cc b/deal.II/source/numerics/point_value_history.cc index a1ea25af97..2cfa0ea914 100644 --- a/deal.II/source/numerics/point_value_history.cc +++ b/deal.II/source/numerics/point_value_history.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2009, 2010, 2012 by Michael Rapson and the deal.II authors +// Copyright (C) 2009, 2010 by Michael Rapson and the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -20,6 +20,8 @@ #include #include +#include + #include @@ -33,10 +35,12 @@ namespace internal /// Only a constructor needed for this class (a struct really) template PointGeometryData - ::PointGeometryData (const std::vector > &new_locations, - const std::vector &new_sol_indices) + ::PointGeometryData (const Point &new_requested_location, + const std::vector > &new_locations, + const std::vector &new_sol_indices) { - locations = new_locations; + requested_location = new_requested_location; + support_point_locations = new_locations; solution_indices = new_sol_indices; } } @@ -47,15 +51,15 @@ namespace internal template PointValueHistory ::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 (); // initialize the std::vector - // make a vector of independent values + // make a vector of independent values independent_values = std::vector > (n_indep, std::vector (0)); } @@ -64,17 +68,17 @@ PointValueHistory template PointValueHistory::PointValueHistory (const DoFHandler & 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 (); // initialize the std::vector - // make a vector for the independent values + // make a vector for the independent values independent_values = std::vector > (n_indep, std::vector (0)); } @@ -87,6 +91,7 @@ PointValueHistory::PointValueHistory (const PointValueHistory & point_value dataset_key = point_value_history.dataset_key; independent_values = point_value_history.independent_values; data_store = point_value_history.data_store; + scalar_field = point_value_history.scalar_field; point_geometry_data = point_value_history.point_geometry_data; pair_data = point_value_history.pair_data; closed = point_value_history.closed; @@ -107,6 +112,7 @@ PointValueHistory::operator= (const PointValueHistory & point_value_history dataset_key = point_value_history.dataset_key; independent_values = point_value_history.independent_values; data_store = point_value_history.data_store; + scalar_field = point_value_history.scalar_field; point_geometry_data = point_value_history.point_geometry_data; pair_data = point_value_history.pair_data; closed = point_value_history.closed; @@ -133,49 +139,49 @@ template void PointValueHistory ::add_point (const Point & 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 > 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 support_point_quadrature (dof_handler->get_fe ().get_unit_support_points ()); FEValues 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::active_cell_iterator cell = dof_handler->begin_active (); typename DoFHandler::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::active_cell_iterator current_cell = cell; std::vector current_fe_index (n_components, 0); // need one index per component fe_values.reinit (cell); @@ -183,32 +189,32 @@ void PointValueHistory 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 test_point - = fe_values.quadrature_point (support_point); + = fe_values.quadrature_point (support_point); if (location.distance (test_point) < - location.distance (current_points [component])) + location.distance (current_points [component])) { - // save the data + // save the data current_points [component] = test_point; current_cell = cell; current_fe_index [component] = support_point; @@ -221,44 +227,45 @@ void PointValueHistory local_dof_indices (dof_handler->get_fe ().dofs_per_cell); std::vector 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 - 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 > >::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 (0)); } @@ -271,46 +278,46 @@ template void PointValueHistory ::add_points (const std::vector > & 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 > 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 support_point_quadrature (dof_handler->get_fe ().get_unit_support_points ()); FEValues 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::active_cell_iterator cell = dof_handler->begin_active (); typename DoFHandler::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 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 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 ::active_cell_iterator > current_cell (locations.size (), cell); fe_values.reinit (cell); @@ -318,8 +325,8 @@ void PointValueHistory std::vector temp_fe_index (n_components, 0); for (unsigned int support_point = 0; support_point < n_support_points; support_point++) { - // setup valid data in the empty - // vectors + // setup valid data in the empty + // vectors unsigned int component = dof_handler->get_fe ().system_to_component_index (support_point).first; temp_points [component] = fe_values.quadrature_point (support_point); temp_fe_index [component] = support_point; @@ -327,8 +334,8 @@ void PointValueHistory std::vector > > current_points (locations.size (), temp_points); // give a valid start point std::vector > current_fe_index (locations.size (), temp_fe_index); - // check each cell to find suitable support - // points + // check each cell to find suitable support + // points for (; cell != endc; cell++) { fe_values.reinit (cell); @@ -341,7 +348,7 @@ void PointValueHistory { if (locations[point].distance (test_point) < locations[point].distance (current_points[point][component])) { - // save the data + // save the data current_points[point][component] = test_point; current_cell[point] = cell; current_fe_index[point][component] = support_point; @@ -361,17 +368,18 @@ void PointValueHistory new_solution_indices.push_back (local_dof_indices[current_fe_index[point][component]]); } - internal::PointValueHistory::PointGeometryData new_point_geometry_data (current_points[point], new_solution_indices); + internal::PointValueHistory::PointGeometryData new_point_geometry_data (locations[point], current_points[point], new_solution_indices); point_geometry_data.push_back (new_point_geometry_data); std::map > >::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 (0)); } @@ -388,22 +396,53 @@ template void PointValueHistory ::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 > vector_size (n_datastreams, std::vector (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 (vector_name, true)); + else + scalar_field.insert (std::pair (vector_name, false)); +} + + + +template +void PointValueHistory +::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 > vector_size (n_datastreams, + std::vector (0)); + pair_data.second = vector_size; + data_store.insert (pair_data); + scalar_field.insert (std::pair (vector_name, true)); // scalar by definition in this method } @@ -440,10 +479,10 @@ void PointValueHistory template template void PointValueHistory -::evaluate_field (const std::string &vector_name, const VECTOR & solution) +::evaluate_field (const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component) { - // must be closed to add data to internal - // members. + // must be closed to add data to internal + // members. Assert (closed, ExcInvalidState ()); Assert (!cleared, ExcInvalidState ()); Assert (n_dofs != 0, ExcDoFHandlerRequired ()); @@ -452,28 +491,297 @@ void PointValueHistory { Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ()); } - typename std::vector >::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 > >::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 >::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 +template +void PointValueHistory +::evaluate_field(const VECTOR & solution, const DataPostprocessor< dim> & data_postprocessor, const Quadrature & 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 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 field_names = data_postprocessor.get_names (); + +// Loop over points and find correct cell + typename std::vector >::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 requested_location = point->requested_location; + typename DoFHandler::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1(), *dof_handler, requested_location).first; + + + fe_values.reinit (cell); + std::vector< Vector< double > > computed_quantities (1, Vector (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 > dummy_normals (1, Point ()); + std::vector > 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 > (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 (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 > dummy_normals (1, Point ()); + std::vector > 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 > (1, evaluation_points[selected_point]), + computed_quantities); + } + + // we now have the data and need to save it +// std::vector component_interpretation = data_postprocessor.get_data_component_interpretation (); + + // loop over data names + typename std::vector::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 > >::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 +template +void PointValueHistory +::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 > >::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 >::iterator point = point_geometry_data.begin (); + Vector 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 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 void PointValueHistory ::start_new_dataset (double key) { - // must be closed to add data to internal - // members. + // must be closed to add data to internal + // members. Assert (closed, ExcInvalidState ()); Assert (!cleared, ExcInvalidState ()); Assert (deep_check (false), ExcDataLostSync ()); @@ -487,8 +795,8 @@ template void PointValueHistory ::push_back_independent (const std::vector &indep_values) { - // must be closed to add data to internal - // members. + // must be closed to add data to internal + // members. Assert (closed, ExcInvalidState ()); Assert (!cleared, ExcInvalidState ()); Assert (indep_values.size () == n_indep, ExcDimensionMismatch (indep_values.size (), n_indep)); @@ -503,13 +811,13 @@ void PointValueHistory template void PointValueHistory -::write_gnuplot (const std::string &base_name) +::write_gnuplot (const std::string &base_name, const std::vector > postprocessor_locations) { AssertThrow (closed, ExcInvalidState ()); AssertThrow (!cleared, ExcInvalidState ()); AssertThrow (deep_check (true), ExcDataLostSync ()); - // write inputs to a file + // write inputs to a file if (n_indep != 0) { std::string filename = base_name + "_indep.gpl"; @@ -517,7 +825,7 @@ void PointValueHistory to_gnuplot << "# Data independent of mesh location\n"; - // write column headings + // write column headings to_gnuplot << "# "; for (unsigned int component = 0; component < n_indep; component++) @@ -526,7 +834,7 @@ void PointValueHistory } to_gnuplot << "\n"; - // write general data stored + // write general data stored for (unsigned int key = 0; key < dataset_key.size (); key++) { to_gnuplot << dataset_key[key]; @@ -543,39 +851,57 @@ void PointValueHistory - // 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 >::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 >::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 > >::iterator - data_store_begin = data_store.begin (); + data_store_begin = data_store.begin (); to_gnuplot << "# "; for (unsigned int component = 0; component < n_indep; component++) @@ -585,18 +911,21 @@ void PointValueHistory 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 > >::iterator - data_store_begin = data_store.begin (); + data_store_begin = data_store.begin (); to_gnuplot << dataset_key[key]; for (unsigned int component = 0; component < n_indep; component++) @@ -606,10 +935,18 @@ void PointValueHistory for (; data_store_begin != data_store.end (); data_store_begin++) { - for (unsigned int component = 0; component < n_components; component++) - { - to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key]; - } + bool scalar = (scalar_field.find (data_store_begin->first))->second; + if (scalar) + { + to_gnuplot << " " << (data_store_begin->second)[data_store_index][key]; + } + else + { + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key]; + } + } } to_gnuplot << "\n"; } @@ -623,48 +960,102 @@ void PointValueHistory template Vector PointValueHistory -::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 dof_vector (dof_handler->n_dofs ()); - typename std::vector >::iterator node = point_geometry_data.begin (); - for (; node != point_geometry_data.end (); node++) + typename std::vector >::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 +Vector PointValueHistory +::mark_locations () +{ + return mark_support_locations (); +} + template void PointValueHistory -::get_points (std::vector > > & locations) +::get_support_locations (std::vector > > & 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 > > actual_points; - typename std::vector >::iterator node = point_geometry_data.begin (); + typename std::vector >::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 +void PointValueHistory +::get_points (std::vector > > & locations) +{ + get_support_locations (locations); +} + + +template +void PointValueHistory +::get_postprocessor_locations (std::vector > & locations, const Quadrature & quadrature) +{ + Assert (!cleared, ExcInvalidState ()); + Assert (n_dofs != 0, ExcDoFHandlerRequired ()); + + locations = std::vector > (); + + FEValues fe_values (dof_handler->get_fe (), quadrature, update_quadrature_points); + unsigned int n_quadrature_points = quadrature.size(); + std::vector > evaluation_points; + +// Loop over points and find correct cell + typename std::vector >::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 requested_location = point->requested_location; + typename DoFHandler::active_cell_iterator cell = GridTools::find_active_cell_around_point (MappingQ1(), *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 void PointValueHistory @@ -675,29 +1066,30 @@ void PointValueHistory out << "Cleared: " << cleared << "\n"; out << "Geometric Data" << "\n"; - typename std::vector >::iterator node = point_geometry_data.begin (); - if (node == point_geometry_data.end ()) + typename std::vector >::iterator point = point_geometry_data.begin (); + if (point == point_geometry_data.end ()) { - out << "No nodes stored currently\n"; + out << "No points stored currently\n"; } else { if (!cleared) - { - out << "# DoF_index : Location (for each component)\n"; - for (; node != point_geometry_data.end (); node++) - { - for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) - { - out << node->solution_indices[component] << " : " << node->locations [component] << "\n"; - } - out << "\n"; - } - } + { + for (; point != point_geometry_data.end (); point++) + { + out << "# Requested location: " << point->requested_location << "\n"; + out << "# DoF_index : Support location (for each component)\n"; + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + out << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n"; + } + out << "\n"; + } + } else - { + { out << "#Cannot access DoF_indices once cleared\n"; - } + } } out << "\n"; @@ -733,8 +1125,8 @@ template bool PointValueHistory ::deep_check (bool strict) { - // test ways that it can fail, if control - // reaches last statement return true + // test ways that it can fail, if control + // reaches last statement return true if (strict) { if (n_indep != 0) @@ -745,19 +1137,19 @@ bool PointValueHistory } } std::map > >::iterator - data_store_begin = data_store.begin (); + data_store_begin = data_store.begin (); if (n_dofs != 0) { for (; data_store_begin != data_store.end (); data_store_begin++) { if ((data_store_begin->second)[0].size () != dataset_key.size ()) return false; - // this loop only tests one - // member for each name, - // i.e. checks the user it will - // not catch internal errors - // which do not update all - // fields for a name. + // this loop only tests one + // member for each name, + // i.e. checks the user it will + // not catch internal errors + // which do not update all + // fields for a name. } } return true; @@ -773,16 +1165,16 @@ bool PointValueHistory if (n_dofs != 0) { std::map > >::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; diff --git a/deal.II/source/numerics/point_value_history.inst.in b/deal.II/source/numerics/point_value_history.inst.in index 9d3fdbeaf1..571307394b 100644 --- a/deal.II/source/numerics/point_value_history.inst.in +++ b/deal.II/source/numerics/point_value_history.inst.in @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2009, 2012 by Michael Rapson and the deal.II authors +// Copyright (C) 2009 by Michael Rapson and the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -19,9 +19,29 @@ for (deal_II_dimension : DIMENSIONS) for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS) -{ +{ template - void PointValueHistory::evaluate_field + void PointValueHistory::evaluate_field (const std::string &, - const VEC &); + const VEC &, + const unsigned int); +} + + +for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS) +{ + template + void PointValueHistory::evaluate_field_at_requested_location + (const std::string &, + const VEC &, + const unsigned int); +} + +for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS) +{ + template + void PointValueHistory::evaluate_field + (const VEC &, + const DataPostprocessor &, + const Quadrature &); } -- 2.39.5