From 561167c68e5a8020f56d18f57c96922468315dda Mon Sep 17 00:00:00 2001 From: rapson Date: Thu, 12 Jul 2012 21:57:18 +0000 Subject: [PATCH] Modified how the PointValueHistory class supports vector fields to allow a more consistent interface. git-svn-id: https://svn.dealii.org/trunk@25695 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/numerics/point_value_history.h | 262 ++-- .../source/numerics/point_value_history.cc | 1125 +++++++++-------- .../numerics/point_value_history.inst.in | 17 +- 3 files changed, 791 insertions(+), 613 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 5d3b023cfa..fc90377f7c 100644 --- a/deal.II/include/deal.II/numerics/point_value_history.h +++ b/deal.II/include/deal.II/numerics/point_value_history.h @@ -101,9 +101,10 @@ namespace internal * just once at the start. However, this method is not suitable for * FiniteElement objects that do not assign dofs to actual mesh locations * (i.e. FEs without @ref GlossSupport "support points") or if adaptive mesh refinement is used. - * The class will throw an exception if the finite element does not have support points - * or any change in the @p dof_handler is identified. (The test is simply - * whether the number of dofs changes.) + * The class will throw an exception if any change to the triangulation is made. Although the + * listeners that detect this change could be used to recompute the nearest support points, + * the location of the support point will most likely change slightly, making the interpretation + * of the data difficult, hence this is not implemented currently. * *
  • Secondly, @p evaluate_field_at_requested_location calls @p * VectorTools::point_value to compute values at the specific point requested. @@ -119,16 +120,21 @@ namespace internal * code using adaptive refinement, but as the output will be meaningful (in * the sense that the quadrature point selected is guaranteed to remain in the * same vicinity, the class does not prevent the use of this method in adaptive - * codes. + * codes. The class provides warnings in the output files if the mesh has changed. * * - * When recording a new mnemonic name, the user must chose whether it should be - * a strictly scalar field (@p add_scalar_field_name) or whether it should - * depend on the @p FE (@p add_field_name). Calling @p add_field_name with a - * @p dof_handler using a vector @p FE will result in a vector field. The vector - * field will be of dimention @p FE.n_components(). Therefore, in this class the - * distinction between vector and scalar is as related to the @p FE and not - * whether a quantity is a logical vector. + * When recording a new mnemonic name, the user must supply a std::vector + * components that lists the (generally FE) components to be extracted from the + * given input. If the user simply wants to extract all the components, or the first + * n components then default parameters and overloaded functions mean that + * the std::vector need not be explicitly supplied to the @p add_field_name method. + * If the @p evaluate_field with a @p DataPostprocessor object is used, the components + * are interpreted as the component of the @p DataPostprocessor return vector. This + * vector can be larger than the FE space. + * + * The class automatically generates names for the data stored based on the mnemoics + * supplied. The methods @p add_component_names and @p add_independent_names allow + * the user to provide lists of names to use instead if desired. * * Following is a little code snippet that shows a common usage of this class: * @@ -184,24 +190,20 @@ class PointValueHistory /** * Constructor linking the class to a - * specific @p DoFHandler. Since this + * specific @p DoFHandler. This * class reads specific data from the @p * DoFHandler and stores it internally * for quick access (in particular dof * indices of closest neighbors to * requested points) the class is fairly * intolerant to changes to the @p - * DoFHandler. Mesh refinement and @p + * DoFHandler if data at support points + * is required. Mesh refinement and @p * DoFRenumbering methods should be * performed before the @p add_points * method is called and adaptive grid - * refinement is not supported (unless it - * is completed before points are added). - * Changes to the @p DoFHandler are - * tested for by looking for a change in - * the number of dofs, which may not - * always occur so the user must be aware - * of these limitations. + * refinement is only supported by some + * methods. * * The user can store extra variables * which do not relate to mesh location @@ -217,9 +219,7 @@ class PointValueHistory const unsigned int n_independent_variables = 0); /** - * Copy constructor explicitly provided, - * although default copy operator would - * be sufficient. This constructor can + * Copy constructor. This constructor can * be safely called with a @p * PointValueHistory object that contains * data, but this could be expensive and @@ -228,9 +228,7 @@ class PointValueHistory PointValueHistory (const PointValueHistory & point_value_history); /** - * Assignment operator explicitly - * provided, although default assignment - * operator would be sufficient. This + * Assignment operator. This * assignment operator can be safely * called once the class is closed and * data added, but this is provided @@ -252,7 +250,7 @@ class PointValueHistory * Add a single point to the class. The * support points (one per component) in * the mesh that are closest to that - * point is found and its details stored + * point are found and their details stored * for use when @p evaluate_field is * called. If more than one point is * required rather used the @p add_points @@ -265,7 +263,7 @@ class PointValueHistory * Add multiple points to the class. The * support points (one per component) in * the mesh that are closest to that - * point is found and its details stored + * point is found and their details stored * for use when @p evaluate_field is * called. If more than one point is * required, rather call this method as it @@ -286,41 +284,59 @@ class PointValueHistory /** * Put another mnemonic string (and hence * @p VECTOR) into the class. This method - * adds storage space for FE.n_components() - * variables and should be used with methods - * that extract all FE components from a vector. + * adds storage space for components.size() + * variables. * This also adds extra entries for points * that are already in the class, so * @p add_field_name and @p add_points can * be called in any order. */ - void add_field_name(const std::string &vector_name); - + void add_field_name(const std::string &vector_name, + const std::vector &components = std::vector ()); + /** * Put another mnemonic string (and hence * @p VECTOR) into the class. This method - * only adds storage space for a single - * variable whether the FE is scalar or not - * so it can be used with - * @p DataPostprocessor objects that compute - * a scalar value(s) from a vector field. + * adds storage space for components.size() + * variables. * This also adds extra entries for points * that are already in the class, so * @p add_field_name and @p add_points can * be called in any order. - */ - void add_scalar_field_name(const std::string &vector_name); + * This method generates a std::vector + * 0, ..., n_components-1 and calls the + * previous function. + */ + void add_field_name(const std::string &vector_name, + const unsigned int n_components); + + /** + * Provide optional names for each component + * of a field. These names will be used + * instead of names generated from the + * field name, if supplied. + */ + void add_component_names(const std::string &vector_name, + const std::vector &component_names); + + /** + * Provide optional names for the + * independent values. These names will + * be used instead of "Indep_...", if + * supplied. + */ + void add_independent_names(const std::vector &independent_names); + + /** * Extract values at the stored points * from the VECTOR supplied and add them * to the new dataset in vector_name. - * The scalar_component input is only - * used if the FE has multiple components - * and the field to be evaluated is scalar. - * In this case, the scalar_component input - * is used to select the component used. If - * a @p DoFHandler is used, one (and only + * The list of components supplied when the + * field was added is used to select + * components to extract. + * If a @p DoFHandler is used, one (and only * one) evaluate_field method * must be called for each dataset (time * step, iteration, etc) for each @@ -328,23 +344,30 @@ class PointValueHistory * ExcDataLostSync error can occur. */ template - void evaluate_field(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0); + void evaluate_field(const std::string &vector_name, const VECTOR & solution); /** - * Compute values using a DataPostprocessor object - * with the VECTOR supplied and add them + * Compute values using a @p DataPostprocessor object + * with the @p VECTOR supplied and add them * to the new dataset in vector_name. - * The DataPostprocessor component_interpretation - * is used to determine select scalar or vector data. - * If the component_is_part_of_vector flag is used, then - * the indicated vector must have the same number of - * components as the FE associated with the dof_handler. - * The DataPostprocessor object is queried for field names. - * In the case of vector fields, the name of the first - * component is used as the overall name of the vector. + * The list of components supplied when the + * field was added is used to select + * components to extract from the @p DataPostprocessor + * return vector. + * This method takes a vector of field names to + * process and is preferred if many fields use the + * same @p DataPostprocessor object as each cell + * is only located once. * The quadrature object supplied is used for all * components of a vector field. + * Although this method will not throw an exception + * if the mesh has changed. (No internal data structures + * are invalidated as the quadrature points are repicked + * each time the function is called.) Nevertheless + * the user must be aware that if the mesh changes the + * point selected will also vary slightly, making + * interpretation of the data more difficult. * If a @p DoFHandler is used, one (and only * one) evaluate_field method * must be called for each dataset (time @@ -353,7 +376,17 @@ class PointValueHistory * ExcDataLostSync error can occur. */ template - void evaluate_field(const VECTOR & solution, const DataPostprocessor & data_postprocessor, const Quadrature & quadrature); + void evaluate_field(const std::vector &vector_names, const VECTOR & solution, const DataPostprocessor & data_postprocessor, const Quadrature & quadrature); + + /** + * Construct a std::vector + * containing only vector_name and + * call the above function. The above function + * is more efficient if multiple fields + * use the same @p DataPostprocessor object. + */ + template + void evaluate_field(const std::string &vector_name, const VECTOR & solution, const DataPostprocessor & data_postprocessor, const Quadrature & quadrature); /** @@ -365,15 +398,11 @@ class PointValueHistory * has been modified because it uses calls to * @p VectorTools::point_value to extract there * data. Therefore, if only this method is used, - * the class can be used with adaptive refinement. - * Specifically, the vector form of - * @p VectorTools::point_value that assumes a - * Q1 mapping is used. - * The scalar_component input is only - * used if the FE has multiple components - * and the field to be evaluated is scalar. - * In this case, the scalar_component input - * is used to select the component used. If + * the class is fully compatible with + * adaptive refinement. + * The list of components supplied when the + * field was added is used to select + * components to extract. If * a @p DoFHandler is used, one (and only * one) evaluate_field method * must be called for each dataset (time @@ -382,7 +411,8 @@ class PointValueHistory * ExcDataLostSync error can occur. */ template - void evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component = 0); + void evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution); // changed + /** * Add the key for the current dataset to @@ -423,6 +453,10 @@ class PointValueHistory * key and independent data. The file * name agrees with the order the points * were added to the class. + * The names of the data columns can + * be supplied using the functions + * @p add_component_names and + * @p add_independent_names. * The support point information is only * meaningful if the dof_handler has not * been changed. Therefore, if adaptive @@ -527,7 +561,7 @@ class PointValueHistory * to the correct number of points by the * method. */ - void get_postprocessor_locations (std::vector > & locations, const Quadrature & quadrature); + void get_postprocessor_locations (const Quadrature & quadrature, std::vector > & locations); /** * Once datasets have been added to the @@ -593,6 +627,19 @@ class PointValueHistory bool deep_check (bool strict); + + /** + * A function that will be triggered + * through signals whenever the + * triangulation is modified. + * + * It is currently used to check + * if the triangulation has changed, + * invalidating precomputed values. + */ + void tria_change_listener (); + + /** * A call has been made to @p * push_back_independent when no @@ -620,26 +667,21 @@ class PointValueHistory /** * A method which requires access to a @p * DoFHandler to be meaningful has been - * called when @p n_dofs is reported to - * be zero (most likely due to default - * constructor being called). Only + * called when @p have_dof_handler is + * false (most likely due to default + * constructor being called). Only * independent variables may be logged * with no DoFHandler. */ DeclException0(ExcDoFHandlerRequired); /** - * @p n_dofs for the @p DoFHandler has - * changed since the constructor was - * called. Possibly the mesh has been - * refined in some way. This suggests + * The triangulation indicated that mesh + * has been refined in some way. This suggests * that the internal dof indices stored * are no longer meaningful. */ - DeclException2(ExcDoFHandlerChanged, - int, - int, - << "Original n_dofs (" << arg1 << ") != current n_dofs (" << arg2 << ")"); + DeclException0(ExcDoFHandlerChanged); private: /** @@ -654,6 +696,14 @@ class PointValueHistory * location. */ std::vector > independent_values; + + /** + * Saves a vector listing component + * names associated with a independent_values. + * This will be an empty vector + * if the user does not supplies names. + */ + std::vector indep_names; /** * Saves data for each mnemonic entry. @@ -667,14 +717,22 @@ class PointValueHistory * FE.n_components () long. */ std::map > > data_store; - + + /** - * Saves a boolean representing whether - * the mnemonic links to a scalar field - * (true) or not (false). + * Saves a vector listing components + * associated with a mnemonic. */ - std::map scalar_field; - + std::map > field_components; + + /** + * Saves a vector listing component + * names associated with a mnemonic. + * This will be an empty vector + * if the user does not supplies names. + */ + std::map > component_names_map; + /** * Saves the location and other mesh info * about support points. @@ -682,11 +740,6 @@ class PointValueHistory std::vector > point_geometry_data; - /** - * A tempory pair used when adding to - * @p data_store. - */ - std::pair > > pair_data; // could possibly be removed /** * Used to enforce @p closed state for some @@ -708,17 +761,26 @@ class PointValueHistory */ 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 - * triangulation has not been refined in - * anyway since the class was set up. - * Refinement will invalidate stored dof - * indices linked to support points. + * Variable to check if the triangulation + * has changed. If it has changed, certain + * data is out of date (especially the + * PointGeometryData::solution_indices. + */ + bool triangulation_changed; + + /** + * A boolean to record whether the class was + * initialized with a DoFHandler or not. */ - unsigned int n_dofs; + bool have_dof_handler; + + /** + * Used to detect signals from the Triangulation. + */ + boost::signals2::connection tria_listener; + /** * Stores the number of independent * variables requested. diff --git a/deal.II/source/numerics/point_value_history.cc b/deal.II/source/numerics/point_value_history.cc index 2cfa0ea914..208cdaabe9 100644 --- a/deal.II/source/numerics/point_value_history.cc +++ b/deal.II/source/numerics/point_value_history.cc @@ -36,10 +36,10 @@ namespace internal template PointGeometryData ::PointGeometryData (const Point &new_requested_location, - const std::vector > &new_locations, - const std::vector &new_sol_indices) + const std::vector > &new_locations, + const std::vector &new_sol_indices) { - requested_location = new_requested_location; + requested_location = new_requested_location; support_point_locations = new_locations; solution_indices = new_sol_indices; } @@ -51,36 +51,45 @@ namespace internal template PointValueHistory ::PointValueHistory (const unsigned int n_independent_variables) : - n_dofs (0), - n_indep (n_independent_variables) + n_indep (n_independent_variables) { closed = false; cleared = false; - // make a vector for keys + triangulation_changed = false; + have_dof_handler = false; + + // make a vector for keys dataset_key = std::vector (); // 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)); + indep_names = std::vector (); } 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_indep (n_independent_variables) { closed = false; cleared = false; - // make a vector to store keys + triangulation_changed = false; + have_dof_handler = true; + + // make a vector to store keys dataset_key = std::vector (); // 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)); + indep_names = std::vector (); + + tria_listener = dof_handler.get_tria().signals.any_change.connect (std_cxx1x::bind (&PointValueHistory::tria_change_listener, + std_cxx1x::ref(*this))); } @@ -90,17 +99,28 @@ PointValueHistory::PointValueHistory (const PointValueHistory & point_value { dataset_key = point_value_history.dataset_key; independent_values = point_value_history.independent_values; + indep_names = point_value_history.indep_names; data_store = point_value_history.data_store; - scalar_field = point_value_history.scalar_field; + field_components = point_value_history.field_components; + component_names_map = point_value_history.component_names_map; point_geometry_data = point_value_history.point_geometry_data; - pair_data = point_value_history.pair_data; + closed = point_value_history.closed; cleared = point_value_history.cleared; dof_handler = point_value_history.dof_handler; - n_dofs = point_value_history.n_dofs; + triangulation_changed = point_value_history.triangulation_changed; + have_dof_handler = point_value_history.have_dof_handler; n_indep = point_value_history.n_indep; + + // What to do with tria_listener? + // Presume subscribe new instance? + if (have_dof_handler) + { + tria_listener = dof_handler->get_tria().signals.any_change.connect (std_cxx1x::bind (&PointValueHistory::tria_change_listener, + std_cxx1x::ref(*this))); + } } @@ -111,17 +131,29 @@ PointValueHistory::operator= (const PointValueHistory & point_value_history { dataset_key = point_value_history.dataset_key; independent_values = point_value_history.independent_values; + indep_names = point_value_history.indep_names; data_store = point_value_history.data_store; - scalar_field = point_value_history.scalar_field; + field_components = point_value_history.field_components; + component_names_map = point_value_history.component_names_map; point_geometry_data = point_value_history.point_geometry_data; - pair_data = point_value_history.pair_data; + closed = point_value_history.closed; cleared = point_value_history.cleared; dof_handler = point_value_history.dof_handler; - n_dofs = point_value_history.n_dofs; + triangulation_changed = point_value_history.triangulation_changed; + have_dof_handler = point_value_history.have_dof_handler; n_indep = point_value_history.n_indep; + + // What to do with tria_listener? + // Presume subscribe new instance? + if (have_dof_handler) + { + tria_listener = dof_handler->get_tria().signals.any_change.connect (std_cxx1x::bind (&PointValueHistory::tria_change_listener, + std_cxx1x::ref(*this))); + } + return * this; } @@ -131,6 +163,10 @@ template PointValueHistory ::~PointValueHistory () { + if (have_dof_handler) + { + tria_listener.disconnect (); + } } @@ -139,49 +175,48 @@ 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 ())); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); - // Implementation assumes that support - // points locations are dofs locations + // Implementation assumes that support + // points locations are dofs locations AssertThrow (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ()); - // FEValues object to extract quadrature - // points from + // FEValues object to extract quadrature + // points from std::vector > 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); @@ -189,32 +224,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; @@ -227,30 +262,30 @@ 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 @@ -261,11 +296,10 @@ void PointValueHistory data_store_begin = data_store.begin (); for (; data_store_begin != data_store.end (); data_store_begin++) { - // add an extra row to each vector - // entry - bool scalar = (scalar_field.find (data_store_begin->first))->second; - for (unsigned int component = 0; - component < (scalar ? 1 : dof_handler->get_fe ().n_components ()); component++) + // add an extra row to each vector + // entry + unsigned int n_components = (field_components.find (data_store_begin->first))->second.size(); + for (unsigned int component = 0; component < n_components; component++) { data_store_begin->second.push_back (std::vector (0)); } @@ -278,46 +312,48 @@ 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 + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); + + + // Implementation assumes that support + // points locations are dofs locations AssertThrow (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ()); - // FEValues object to extract quadrature - // points from + // FEValues object to extract quadrature + // points from std::vector > 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); @@ -325,8 +361,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; @@ -334,8 +370,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); @@ -348,7 +384,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; @@ -373,13 +409,13 @@ void PointValueHistory 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 - bool scalar = (scalar_field.find (data_store_begin->first))->second; - for (unsigned int component = 0; component < (scalar ? 1 : dof_handler->get_fe ().n_components ()); component++) + // add an extra row to each vector + // entry + unsigned int n_components = (field_components.find (data_store_begin->first))->second.size(); + for (unsigned int component = 0; component < n_components; component++) { data_store_begin->second.push_back (std::vector (0)); } @@ -394,58 +430,84 @@ void PointValueHistory template void PointValueHistory -::add_field_name (const std::string &vector_name) +::add_field_name (const std::string &vector_name, const std::vector &components) { - // can't be closed to add additional points - // or vectors + // can't be closed to add additional points + // or vectors AssertThrow (!closed, ExcInvalidState ()); - AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); AssertThrow (!cleared, ExcInvalidState ()); - AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); - // make and add a new vector - // point_geometry_data.size() long + // Make a field_components list, if not supplied + std::vector temp_components = components; + if (temp_components.size() == 0) + { + // Generate component list + // 0,..,fe.n_components -1 + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + temp_components.push_back (component); + } + field_components.insert (std::pair > (vector_name, temp_components)); + + // insert an empty vector of strings + // to ensure each field has an entry + // in the map + std::pair > empty_names (vector_name, std::vector ()); + component_names_map.insert (empty_names); + + // make and add a new vector + // point_geometry_data.size() long + std::pair > > pair_data; pair_data.first = vector_name; - int n_datastreams = point_geometry_data.size () * (dof_handler->get_fe ().n_components ()); // each point has n_components sub parts + int n_datastreams = point_geometry_data.size () * temp_components.size(); // each point has field_components sub parts std::vector < std::vector > 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) +::add_field_name(const std::string &vector_name, const unsigned int n_components) { - // can't be closed to add additional points - // or vectors - AssertThrow (!closed, ExcInvalidState ()); - AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); - AssertThrow (!cleared, ExcInvalidState ()); - AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + std::vector temp_components; + // Generate component list + // 0,..,n_components -1 + for (unsigned int component = 0; component < n_components; component++) + temp_components.push_back (component); + + add_field_name (vector_name, temp_components); +} - // make and add a new vector - // point_geometry_data.size() long - pair_data.first = vector_name; - int n_datastreams = point_geometry_data.size (); // each point has 1 value for this function - std::vector < std::vector > 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 +template +void PointValueHistory +::add_component_names(const std::string &vector_name, + const std::vector &component_names) +{ + typename std::map >::iterator names = component_names_map.find(vector_name); + Assert (names != component_names_map.end(), ExcMessage("vector_name not in class")); + + typename std::map >::iterator components = field_components.find(vector_name); + Assert (components != field_components.end(), ExcMessage("vector_name not in class")); + Assert (component_names.size() == components->second.size(), ExcDimensionMismatch (component_names.size(), components->second.size())); + + names->second = component_names; } +template +void PointValueHistory +::add_independent_names(const std::vector &independent_names) +{ + Assert (independent_names.size() == n_indep, ExcDimensionMismatch (independent_names.size(), n_indep)); + + indep_names = independent_names; +} + template void PointValueHistory @@ -462,6 +524,7 @@ void PointValueHistory { cleared = true; dof_handler = 0; + have_dof_handler = false; } // Need to test that the internal data has a full and complete dataset for @@ -479,54 +542,49 @@ void PointValueHistory template template void PointValueHistory -::evaluate_field (const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component) +::evaluate_field (const std::string &vector_name, const VECTOR & solution) { - // must be closed to add data to internal - // members. + // must be closed to add data to internal + // members. Assert (closed, ExcInvalidState ()); Assert (!cleared, ExcInvalidState ()); - Assert (n_dofs != 0, ExcDoFHandlerRequired ()); - Assert (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); + if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0 { Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ()); } - // Look up the field name and get an - // iterator for the map. Doing this - // up front means that it only needs - // to be done once and also allows us - // to check vector_name is in the map. + // Look up the field name and get an + // iterator for the map. Doing this + // up front means that it only needs + // to be done once and also allows us + // to check vector_name is in the map. typename std::map > >::iterator data_store_field = data_store.find(vector_name); Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class")); + // Repeat for field_components + typename std::map >::iterator components = field_components.find(vector_name); + Assert (components != field_components.end(), ExcMessage("vector_name not in class")); + + Assert (components->second.size () <= dof_handler->get_fe ().n_components (), ExcIndexRange (components->second.size (), 0, dof_handler->get_fe ().n_components () + 1)); // 0 components legal as request to generate list internally! + for (unsigned int component = 0; component < components->second.size (); component++) + { + Assert (components->second [component] < dof_handler->get_fe ().n_components (), ExcIndexRange (components->second[component], 0, dof_handler->get_fe ().n_components ())); + } // hopefully the compiler will optimise + typename std::vector >::iterator point = point_geometry_data.begin (); for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++) + { + // Look up the component to add + // in field_components, and + // access the data associated with + // that component + for (unsigned int comp = 0; comp < components->second.size(); comp++) { - if (!(scalar_field.find (vector_name))->second) - { - // The solution_indices data has a - // one-to-one correspondance to FE - // components therefore we simply - // step through each point, to access - // the Solution_index and the vector - // index - for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) - { - unsigned int solution_index = point->solution_indices[component]; - data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (solution (solution_index)); - } - } - else - { - // We have found a scalar field - // added to a dof_handler with a - // vector FE. Therefore, we take - // the value specified by the user - // in the scalar_component input. - Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError()); - unsigned int solution_index = point->solution_indices[scalar_component]; - data_store_field->second[data_store_index].push_back (solution (solution_index)); - } + unsigned int solution_index = point->solution_indices[components->second[comp]]; + data_store_field->second[data_store_index * components->second.size() + comp].push_back (solution (solution_index)); } + } } @@ -536,243 +594,225 @@ void PointValueHistory template template void PointValueHistory -::evaluate_field(const VECTOR & solution, const DataPostprocessor< dim> & data_postprocessor, const Quadrature & quadrature) +::evaluate_field(const std::vector &vector_names, const VECTOR & solution, const DataPostprocessor< dim> & data_postprocessor, const Quadrature & quadrature) { - // must be closed to add data to internal - // members. + // must be closed to add data to internal + // members. Assert (closed, ExcInvalidState ()); Assert (!cleared, ExcInvalidState ()); - Assert (n_dofs != 0, ExcDoFHandlerRequired ()); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0 { Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ()); } -// Make an FEValues object + // Make an FEValues object const UpdateFlags update_flags = data_postprocessor.get_needed_update_flags() | update_quadrature_points; Assert (!(update_flags & update_normal_vectors), - ExcMessage("The update of normal vectors may not be requested for evaluation of " - "data on cells via DataPostprocessor.")); + ExcMessage("The update of normal vectors may not be requested for evaluation of " + "data on cells via DataPostprocessor.")); FEValues 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 (); + unsigned int n_output_variables = data_postprocessor.get_names().size(); -// Loop over points and find correct cell + // Loop over points and find correct cell typename std::vector >::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 ++; - } - } + // 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 (n_output_variables)); // just one point needed + + // The case of a scalar FE + if (n_components == 1) + { + // Extract data for the + // PostProcessor object + std::vector< double > uh (n_quadrature_points, 0.0); + std::vector< Tensor< 1, dim > > duh (n_quadrature_points, Tensor <1, dim> ()); + std::vector< Tensor< 2, dim > > dduh (n_quadrature_points, Tensor <2, dim> ()); + std::vector > 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 + // loop over data names + typename std::vector::const_iterator name = vector_names.begin(); + for (; name != vector_names.end(); name++) + { + typename std::map > >::iterator data_store_field = data_store.find(*name); + Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class")); + // Repeat for field_components + typename std::map >::iterator components = field_components.find(*name); + Assert (components != field_components.end(), ExcMessage("vector_name not in class")); + + // Push back components.size() computed + // quantities. Each component requested + // must be less than the postprocessor + // n_output_variables. + for (unsigned int comp = 0; comp < components->second.size(); comp++) + { + Assert (components->second[comp] < n_output_variables, ExcIndexRange (components->second[comp], 0, n_output_variables)); + data_store_field->second[data_store_index * components->second.size() + comp].push_back (computed_quantities[0](components->second[comp])); + } + } } // end of loop over points } - +template +template +void PointValueHistory +::evaluate_field(const std::string &vector_name, const VECTOR & solution, const DataPostprocessor & data_postprocessor, const Quadrature & quadrature) +{ + std::vector vector_names; + vector_names.push_back (vector_name); + evaluate_field (vector_names, solution, data_postprocessor, quadrature); +} template template void PointValueHistory -::evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution, const unsigned int scalar_component) +::evaluate_field_at_requested_location(const std::string &vector_name, const VECTOR & solution) { - // must be closed to add data to internal - // members. + // must be closed to add data to internal + // members. Assert (closed, ExcInvalidState ()); Assert (!cleared, ExcInvalidState ()); - Assert (n_dofs != 0, ExcDoFHandlerRequired ()); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0 { Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ()); } - // Look up the field name and get an - // iterator for the map. Doing this - // up front means that it only needs - // to be done once and also allows us - // to check vector_name is in the map. + // Look up the field name and get an + // iterator for the map. Doing this + // up front means that it only needs + // to be done once and also allows us + // to check vector_name is in the map. typename std::map > >::iterator data_store_field = data_store.find(vector_name); Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class")); + // Repeat for field_components + typename std::map >::iterator components = field_components.find(vector_name); + Assert (components != field_components.end(), ExcMessage("vector_name not in class")); + + Assert (components->second.size () <= dof_handler->get_fe ().n_components (), ExcIndexRange (components->second.size (), 0, dof_handler->get_fe ().n_components () + 1)); // 0 components legal as request to generate list internally! + for (unsigned int component = 0; component < components->second.size (); component++) + { + Assert (components->second [component] < dof_handler->get_fe ().n_components (), ExcIndexRange (components->second[component], 0, dof_handler->get_fe ().n_components ())); + } // hopefully the compiler will optimise + + typename std::vector >::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); + + + // Look up the component to add + // in field_components, and + // access the data associated with + // that component + for (unsigned int comp = 0; comp < components->second.size(); comp++) { - // Make a Vector for the value - // at the point. It will have as many - // components as there are in the fe. - VectorTools::point_value (*dof_handler, solution, point->requested_location, value); - - if (!(scalar_field.find (vector_name))->second) - { - // The solution_indices data has a - // one-to-one correspondance to FE - // components therefore we simply - // step through each point, to access - // the Solution_index and the vector - // index - for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) - { - data_store_field->second[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (value (component)); - } - } - else - { - // We have found a scalar field - // added to a dof_handler with a - // vector FE. Therefore, we take - // the value specified by the user - // in the scalar_component input. - Assert (data_store_field->second.size() == point_geometry_data.size (), ExcInternalError()); - data_store_field->second[data_store_index].push_back (value (scalar_component)); - } + data_store_field->second[data_store_index * components->second.size() + comp].push_back (value (components->second[comp])); } + } } @@ -780,8 +820,8 @@ 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 ()); @@ -795,8 +835,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)); @@ -817,7 +857,7 @@ void PointValueHistory 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"; @@ -825,16 +865,26 @@ 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++) - { - to_gnuplot << " "; - } - to_gnuplot << "\n"; - - // write general data stored + if (indep_names.size() > 0) + { + for (unsigned int name = 0; name < indep_names.size(); name++) + { + to_gnuplot << "<" << indep_names [name] << "> "; + } + to_gnuplot << "\n"; + } + else + { + for (unsigned int component = 0; component < n_indep; component++) + { + to_gnuplot << " "; + } + to_gnuplot << "\n"; + } + // write general data stored for (unsigned int key = 0; key < dataset_key.size (); key++) { to_gnuplot << dataset_key[key]; @@ -851,81 +901,104 @@ void PointValueHistory - // write points to a file - if (n_dofs != 0) + // write points to a file + if (have_dof_handler) { - AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); - AssertThrow (postprocessor_locations.size() == 0 || postprocessor_locations.size() == point_geometry_data.size(), ExcDimensionMismatch (postprocessor_locations.size(), point_geometry_data.size())); - // We previously required the - // number of dofs to remain the - // same to provide some sort of - // test on the relevance of the - // support point indices stored. - // We now relax that to allow - // adaptive refinement strategies - // to make use of the - // evaluate_field_requested_locations - // method. Note that the support point - // information is not meaningful if - // the number of dofs has changed. - //AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); - - unsigned int n_components = dof_handler->get_fe ().n_components (); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + AssertThrow (postprocessor_locations.size() == 0 || postprocessor_locations.size() == point_geometry_data.size(), ExcDimensionMismatch (postprocessor_locations.size(), point_geometry_data.size())); + // We previously required the + // number of dofs to remain the + // same to provide some sort of + // test on the relevance of the + // support point indices stored. + // We now relax that to allow + // adaptive refinement strategies + // to make use of the + // evaluate_field_requested_locations + // method. Note that the support point + // information is not meaningful if + // the number of dofs has changed. + //AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); typename std::vector >::iterator point = point_geometry_data.begin (); for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++) { - // for each point, open a file to - // be written to + // for each point, open a file to + // be written to std::string filename = base_name + "_" + Utilities::int_to_string (data_store_index, 2) + ".gpl"; // store by order pushed back - // due to - // Utilities::int_to_string(data_store_index, - // 2) call, can handle up to 100 - // points + // due to + // Utilities::int_to_string(data_store_index, + // 2) call, can handle up to 100 + // points std::ofstream to_gnuplot (filename.c_str ()); - // put helpful info about the - // support point into the file as - // comments - to_gnuplot << "# Requested location: " << point->requested_location << "\n"; + // put helpful info about the + // support point into the file as + // comments + to_gnuplot << "# Requested location: " << point->requested_location << "\n"; to_gnuplot << "# DoF_index : Support location (for each component)\n"; - for (unsigned int component = 0; component < n_components; component++) + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) { to_gnuplot << "# " << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n"; } - if (postprocessor_locations.size() != 0) - { - to_gnuplot << "# Postprocessor location: " << postprocessor_locations[data_store_index] << "\n"; - } + if (triangulation_changed) + to_gnuplot << "# (Original components and locations, may be invalidated by mesh change.)\n"; + + if (postprocessor_locations.size() != 0) + { + to_gnuplot << "# Postprocessor location: " << postprocessor_locations[data_store_index]; + if (triangulation_changed) + to_gnuplot << " (may be approximate)\n"; + } to_gnuplot << "#\n"; - // write column headings + + // write column headings std::map > >::iterator - data_store_begin = data_store.begin (); + data_store_begin = data_store.begin (); to_gnuplot << "# "; - for (unsigned int component = 0; component < n_indep; component++) - { - to_gnuplot << " "; - } + if (indep_names.size() > 0) + { + for (unsigned int name = 0; name < indep_names.size(); name++) + { + to_gnuplot << "<" << indep_names [name] << "> "; + } + } + else + { + for (unsigned int component = 0; component < n_indep; component++) + { + to_gnuplot << " "; + } + } for (; data_store_begin != data_store.end (); data_store_begin++) { - if (!(scalar_field.find (data_store_begin->first))->second) - for (unsigned int component = 0; component < n_components; component++) - { - to_gnuplot << "<" << data_store_begin->first << "_" << component << "> "; - } - else - to_gnuplot << "<" << data_store_begin->first << "_" << 0 << "> "; + unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components + std::vector names = (component_names_map.find (data_store_begin->first))->second; + if (names.size() > 0) + { + for (unsigned int component = 0; component < names.size(); component++) + { + to_gnuplot << "<" << names[component] << "> "; + } + } + else + { + for (unsigned int component = 0; component < n_comps; component++) + { + to_gnuplot << "<" << data_store_begin->first << "_" << component << "> "; + } + } } to_gnuplot << "\n"; - // write data stored for the point + // write data stored for the point for (unsigned int key = 0; key < dataset_key.size (); key++) { std::map > >::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++) @@ -935,18 +1008,12 @@ void PointValueHistory for (; data_store_begin != data_store.end (); data_store_begin++) { - bool scalar = (scalar_field.find (data_store_begin->first))->second; - if (scalar) - { - to_gnuplot << " " << (data_store_begin->second)[data_store_index][key]; - } - else - { - for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) - { - to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key]; - } - } + unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components + + for (unsigned int component = 0; component < n_comps; component++) + { + to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_comps + component][key]; + } } to_gnuplot << "\n"; } @@ -962,11 +1029,11 @@ template Vector PointValueHistory ::mark_support_locations () { - // a method to put a one at each point on - // the grid where a location is defined - AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); + // a method to put a one at each point on + // the grid where a location is defined AssertThrow (!cleared, ExcInvalidState ()); - AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); Vector dof_vector (dof_handler->n_dofs ()); @@ -986,7 +1053,7 @@ template Vector PointValueHistory ::mark_locations () { - return mark_support_locations (); + return mark_support_locations (); } @@ -994,9 +1061,9 @@ template void PointValueHistory ::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 ())); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); + AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); std::vector > > actual_points; typename std::vector >::iterator point = point_geometry_data.begin (); @@ -1013,16 +1080,16 @@ template void PointValueHistory ::get_points (std::vector > > & locations) { - get_support_locations (locations); + get_support_locations (locations); } template void PointValueHistory -::get_postprocessor_locations (std::vector > & locations, const Quadrature & quadrature) +::get_postprocessor_locations (const Quadrature & quadrature, std::vector > & locations) { Assert (!cleared, ExcInvalidState ()); - Assert (n_dofs != 0, ExcDoFHandlerRequired ()); + AssertThrow (have_dof_handler, ExcDoFHandlerRequired ()); locations = std::vector > (); @@ -1030,29 +1097,30 @@ void PointValueHistory unsigned int n_quadrature_points = quadrature.size(); std::vector > evaluation_points; -// Loop over points and find correct cell + // 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]); + // 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]); } } @@ -1064,6 +1132,8 @@ void PointValueHistory out << "***PointValueHistory status output***\n\n"; out << "Closed: " << closed << "\n"; out << "Cleared: " << cleared << "\n"; + out << "Triangulation_changed: " << triangulation_changed << "\n"; + out << "Have_dof_handler: " << have_dof_handler << "\n"; out << "Geometric Data" << "\n"; typename std::vector >::iterator point = point_geometry_data.begin (); @@ -1074,28 +1144,37 @@ void PointValueHistory else { if (!cleared) - { - for (; point != point_geometry_data.end (); point++) - { - out << "# Requested location: " << point->requested_location << "\n"; - out << "# DoF_index : Support location (for each component)\n"; - for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) - { - out << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n"; - } - out << "\n"; - } - } + { + for (; point != point_geometry_data.end (); point++) + { + out << "# Requested location: " << point->requested_location << "\n"; + out << "# DoF_index : Support location (for each component)\n"; + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + out << point->solution_indices[component] << " : " << point->support_point_locations [component] << "\n"; + } + out << "\n"; + } + } else - { + { out << "#Cannot access DoF_indices once cleared\n"; - } + } } out << "\n"; if (independent_values.size () != 0) { out << "Independent value(s): " << independent_values.size () << " : " << independent_values[0].size () << "\n"; + if (indep_names.size() > 0) + { + out << "Names: "; + for (unsigned int name = 0; name < indep_names.size(); name++) + { + out << "<" << indep_names [name] << "> "; + } + out << "\n"; + } } else { @@ -1106,15 +1185,35 @@ void PointValueHistory data_store_begin = data_store.begin (); for (; data_store_begin != data_store.end (); data_store_begin++) { + // Find field mnemonic + std::string vector_name = data_store_begin->first; + typename std::map >::iterator components = field_components.find(vector_name); + Assert (components != field_components.end(), ExcMessage("vector_name not in class")); + typename std::map >::iterator component_names = component_names_map.find(vector_name); + Assert (component_names != component_names_map.end(), ExcMessage("vector_name not in class")); + if (data_store_begin->second.size () != 0) { - out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << (data_store_begin->second)[0].size () << "\n"; + out << data_store_begin->first << ": " << data_store_begin->second.size () << " ("; + out << components->second.size() << ") : "; + out << (data_store_begin->second)[0].size () << "\n"; } else { - out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << "No points added" << "\n"; + out << data_store_begin->first << ": " << data_store_begin->second.size () << " ("; + out << components->second.size() << ") : "; + out << "No points added" << "\n"; } - } + // add names, if available + if (component_names->second.size() > 0) + { + for (unsigned int name = 0; name < component_names->second.size(); name++) + { + out << "<" << component_names->second[name] << "> "; + } + out << "\n"; + } + } out << "\n"; out << "***end of status output***\n\n"; } @@ -1125,8 +1224,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) @@ -1137,19 +1236,19 @@ bool PointValueHistory } } std::map > >::iterator - data_store_begin = data_store.begin (); - if (n_dofs != 0) + data_store_begin = data_store.begin (); + if (have_dof_handler) { for (; data_store_begin != data_store.end (); data_store_begin++) { if ((data_store_begin->second)[0].size () != dataset_key.size ()) return false; - // this loop only tests one - // member for each name, - // i.e. checks the user it will - // not catch internal errors - // which do not update all - // fields for a name. + // this loop only tests one + // member for each name, + // i.e. checks the user it will + // not catch internal errors + // which do not update all + // fields for a name. } } return true; @@ -1162,19 +1261,19 @@ bool PointValueHistory } } - if (n_dofs != 0) + if (have_dof_handler) { 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; @@ -1182,6 +1281,14 @@ bool PointValueHistory +template +void PointValueHistory +::tria_change_listener () +{ + triangulation_changed = true; +} + + // explicit instantiations #include "point_value_history.inst" diff --git a/deal.II/source/numerics/point_value_history.inst.in b/deal.II/source/numerics/point_value_history.inst.in index 571307394b..e175ccbb23 100644 --- a/deal.II/source/numerics/point_value_history.inst.in +++ b/deal.II/source/numerics/point_value_history.inst.in @@ -23,8 +23,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS) template void PointValueHistory::evaluate_field (const std::string &, - const VEC &, - const unsigned int); + const VEC &); } @@ -33,15 +32,25 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS) template void PointValueHistory::evaluate_field_at_requested_location (const std::string &, + const VEC &); +} + +for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS) +{ + template + void PointValueHistory::evaluate_field + (const std::vector &, const VEC &, - const unsigned int); + const DataPostprocessor &, + const Quadrature &); } for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS) { template void PointValueHistory::evaluate_field - (const VEC &, + (const std::string &, + const VEC &, const DataPostprocessor &, const Quadrature &); } -- 2.39.5