From 0ce0dbc021a6580c1ede5ec86e2bb79818d20ffd Mon Sep 17 00:00:00 2001 From: rapson Date: Tue, 7 Aug 2012 21:31:47 +0000 Subject: [PATCH] Changed the method of specifying components from a list to a component mask. Documented above changes, and provided comments clarifying that GridTools::find_active_cell_around_point can be used in the add_point(s) methods. White space changes to improve the code structure. git-svn-id: https://svn.dealii.org/trunk@25763 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/numerics/point_value_history.h | 74 +++--- .../source/numerics/point_value_history.cc | 220 ++++++++++-------- 2 files changed, 161 insertions(+), 133 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 1db7ab4041..73152cd852 100644 --- a/deal.II/include/deal.II/numerics/point_value_history.h +++ b/deal.II/include/deal.II/numerics/point_value_history.h @@ -143,14 +143,17 @@ namespace internal * course at the expense of getting only an approximate result.) * * - * When recording a new mnemonic name, the user must supply a std::vector - * components that lists the @ref GlossComponent "(vector) 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. + * When recording a new mnemonic name, the user must supply a std::vector + * component_mask to indicate the @ref GlossComponent "(vector) components" to be + * extracted from the given input. If the user simply wants to extract all the components, + * then default parameters + * the mask 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 component_mask + * interpreted as the mask of the @p DataPostprocessor return vector. The size of this + * mask can be different to that of the FE space, but must be provided when the + * @p add_field_name method is called. One variant of the @p add_field_name method allows + * provides an unsigned int input to construct a suitable mask, if all values from the + * @p DataPostprocessor are desired. * * The class automatically generates names for the data stored based on the mnemonics * supplied. The methods @p add_component_names and @p add_independent_names allow @@ -304,20 +307,21 @@ class PointValueHistory /** * Put another mnemonic string (and hence * @p VECTOR) into the class. This method - * adds storage space for components.size() - * variables. + * adds storage space for variables equal + * to the number of true values in + * component_mask. * 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, - const std::vector &components = std::vector ()); + const std::vector &component_mask = std::vector ()); /** * Put another mnemonic string (and hence * @p VECTOR) into the class. This method - * adds storage space for components.size() + * adds storage space for n_components * variables. * This also adds extra entries for points * that are already in the class, so @@ -353,7 +357,7 @@ class PointValueHistory * Extract values at the stored points * from the VECTOR supplied and add them * to the new dataset in vector_name. - * The list of components supplied when the + * The component mask supplied when the * field was added is used to select * components to extract. * If a @p DoFHandler is used, one (and only @@ -365,14 +369,14 @@ class PointValueHistory */ template void evaluate_field(const std::string &name, - const VECTOR & solution); + const VECTOR & solution); /** * Compute values using a @p DataPostprocessor object * with the @p VECTOR supplied and add them * to the new dataset in vector_name. - * The list of components supplied when the + * The component_mask supplied when the * field was added is used to select * components to extract from the @p DataPostprocessor * return vector. @@ -398,9 +402,9 @@ class PointValueHistory */ template void evaluate_field(const std::vector &names, - const VECTOR & solution, - const DataPostprocessor & data_postprocessor, - const Quadrature & quadrature); + const VECTOR & solution, + const DataPostprocessor & data_postprocessor, + const Quadrature & quadrature); /** * Construct a std::vector @@ -411,9 +415,9 @@ class PointValueHistory */ template void evaluate_field(const std::string &name, - const VECTOR & solution, - const DataPostprocessor & data_postprocessor, - const Quadrature & quadrature); + const VECTOR & solution, + const DataPostprocessor & data_postprocessor, + const Quadrature & quadrature); /** @@ -427,7 +431,7 @@ class PointValueHistory * data. Therefore, if only this method is used, * the class is fully compatible with * adaptive refinement. - * The list of components supplied when the + * The component_mask supplied when the * field was added is used to select * components to extract. If * a @p DoFHandler is used, one (and only @@ -439,7 +443,7 @@ class PointValueHistory */ template void evaluate_field_at_requested_location(const std::string &name, - const VECTOR & solution); + const VECTOR & solution); /** @@ -502,7 +506,7 @@ class PointValueHistory * locations output. */ void write_gnuplot (const std::string &base_name, - const std::vector > postprocessor_locations = std::vector > ()); + const std::vector > postprocessor_locations = std::vector > ()); /** @@ -545,8 +549,8 @@ class PointValueHistory * The function mark_support_locations replaces * it and reflects the fact that the locations * marked are actually the support points. - * - * @deprecated + * + * @deprecated */ Vector mark_locations(); @@ -593,7 +597,7 @@ class PointValueHistory * method. */ void get_postprocessor_locations (const Quadrature & quadrature, - std::vector > & locations); + std::vector > & locations); /** * Once datasets have been added to the @@ -649,8 +653,8 @@ class PointValueHistory * Check the internal data sizes to test * for a loss of data sync. This is often * used in @p Assert statements with the - * @p ExcDataLostSync exception. If @p - * strict is @p false this method returns + * @p ExcDataLostSync exception. If @p strict + * is @p false this method returns * @p true if all sizes are within 1 of * each other (needed to allow data to be * added), with @p strict = @p true they @@ -695,7 +699,7 @@ class PointValueHistory DeclException0(ExcDoFHandlerRequired); /** - * The triangulation indicated that mesh + * The triangulation indicated that mesh * has been refined in some way. This suggests * that the internal dof indices stored * are no longer meaningful. @@ -737,12 +741,12 @@ class PointValueHistory */ std::map > > data_store; + /** + * Saves a component mask + * for each mnemonic. + */ + std::map > component_mask; - /** - * Saves a vector listing components - * associated with a mnemonic. - */ - std::map > field_components; /** * Saves a vector listing component diff --git a/deal.II/source/numerics/point_value_history.cc b/deal.II/source/numerics/point_value_history.cc index 57fdb4fb28..a708e95cd4 100644 --- a/deal.II/source/numerics/point_value_history.cc +++ b/deal.II/source/numerics/point_value_history.cc @@ -24,6 +24,8 @@ #include +#include + DEAL_II_NAMESPACE_OPEN @@ -101,7 +103,7 @@ PointValueHistory::PointValueHistory (const PointValueHistory & point_value independent_values = point_value_history.independent_values; indep_names = point_value_history.indep_names; data_store = point_value_history.data_store; - field_components = point_value_history.field_components; + component_mask = point_value_history.component_mask; component_names_map = point_value_history.component_names_map; point_geometry_data = point_value_history.point_geometry_data; @@ -133,7 +135,7 @@ PointValueHistory::operator= (const PointValueHistory & point_value_history independent_values = point_value_history.independent_values; indep_names = point_value_history.indep_names; data_store = point_value_history.data_store; - field_components = point_value_history.field_components; + component_mask = point_value_history.component_mask; component_names_map = point_value_history.component_names_map; point_geometry_data = point_value_history.point_geometry_data; @@ -234,6 +236,14 @@ void PointValueHistory // check each cell to find a suitable // support points + // GridTools::find_active_cell_around_point + // is an alternative. That method is not + // used here mostly because of the history + // of the class. The algorithm used in + // add_points below may be slightly more + // efficient than find_active_cell_around_point + // because it operates on a set of points. + for (; cell != endc; cell++) { fe_values.reinit (cell); @@ -276,7 +286,7 @@ void PointValueHistory // will still lie within one cell distance. // calling // GridTools::find_active_cell_around_point - // to obtain a cell to search may be an + // to obtain a cell to search is an // option for these methods, but currently // the GridTools method does not cater for // a vector of points, and does not seem to @@ -298,8 +308,9 @@ void PointValueHistory { // 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++) + std::vector current_mask = (component_mask.find (data_store_begin->first))->second; + unsigned int n_stored = std::count(current_mask.begin(), current_mask.end(), true); + for (unsigned int component = 0; component < n_stored; component++) { data_store_begin->second.push_back (std::vector (0)); } @@ -372,6 +383,13 @@ void PointValueHistory // check each cell to find suitable support // points + // GridTools::find_active_cell_around_point + // is an alternative. That method is not + // used here mostly because of the history + // of the class. The algorithm used here + // may be slightly more + // efficient than find_active_cell_around_point + // because it operates on a set of points. for (; cell != endc; cell++) { fe_values.reinit (cell); @@ -414,8 +432,9 @@ void PointValueHistory { // 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++) + std::vector current_mask = (component_mask.find (data_store_begin->first))->second; + unsigned int n_stored = std::count(current_mask.begin(), current_mask.end(), true); + for (unsigned int component = 0; component < n_stored; component++) { data_store_begin->second.push_back (std::vector (0)); } @@ -430,7 +449,7 @@ void PointValueHistory template void PointValueHistory -::add_field_name (const std::string &vector_name, const std::vector &components) +::add_field_name (const std::string &vector_name, const std::vector &mask) { // can't be closed to add additional points // or vectors @@ -440,16 +459,13 @@ void PointValueHistory AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); - // Make a field_components list, if not supplied - std::vector temp_components = components; - if (temp_components.size() == 0) + // Make a component_mask, if not supplied + std::vector temp_mask = mask; + if (temp_mask.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); + temp_mask = std::vector (dof_handler->get_fe().n_components(), true); } - field_components.insert (std::pair > (vector_name, temp_components)); + component_mask.insert (std::pair > (vector_name, temp_mask)); // insert an empty vector of strings // to ensure each field has an entry @@ -461,7 +477,8 @@ void PointValueHistory // point_geometry_data.size() long std::pair > > pair_data; pair_data.first = vector_name; - int n_datastreams = point_geometry_data.size () * temp_components.size(); // each point has field_components sub parts + unsigned int n_stored = std::count(temp_mask.begin(), temp_mask.end(), true); + int n_datastreams = point_geometry_data.size () * n_stored; // each point has n_stored sub parts std::vector < std::vector > vector_size (n_datastreams, std::vector (0)); pair_data.second = vector_size; @@ -473,13 +490,8 @@ template void PointValueHistory ::add_field_name(const std::string &vector_name, const unsigned int n_components) { - 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); + std::vector temp_mask (n_components, true); + add_field_name (vector_name, temp_mask); } @@ -491,9 +503,10 @@ void PointValueHistory 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())); + typename std::map >::iterator mask = component_mask.find(vector_name); + Assert (mask != component_mask.end(), ExcMessage("vector_name not in class")); + unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true); + Assert (component_names.size() == n_stored, ExcDimensionMismatch (component_names.size(), n_stored)); names->second = component_names; } @@ -562,27 +575,29 @@ void PointValueHistory // 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")); + // Repeat for component_mask + typename std::map >::iterator mask = component_mask.find(vector_name); + Assert (mask != component_mask.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 + Assert (mask->second.size () == dof_handler->get_fe ().n_components (), ExcDimensionMismatch (mask->second.size (), dof_handler->get_fe ().n_components ())); + unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true); 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 + // Look up the components to add + // in the component_mask, and // access the data associated with - // that component - for (unsigned int comp = 0; comp < components->second.size(); comp++) + // those components + + for (unsigned int store_index = 0, comp = 0; comp < mask->second.size(); comp++) { - 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)); + if (mask->second[comp]) + { + unsigned int solution_index = point->solution_indices[comp]; + data_store_field->second[data_store_index * n_stored + store_index].push_back (solution (solution_index)); + store_index++; + } } } } @@ -731,18 +746,22 @@ void PointValueHistory { 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++) + // Repeat for component_mask + typename std::map >::iterator mask = component_mask.find(*name); + Assert (mask != component_mask.end(), ExcMessage("vector_name not in class")); + + Assert (mask->second.size () == n_output_variables, ExcDimensionMismatch (mask->second.size (),n_output_variables)); + unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true); + + // Push back computed quantities according + // to the component_mask. + for (unsigned int store_index = 0, comp = 0; comp < mask->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])); + if (mask->second[comp]) + { + data_store_field->second[data_store_index * n_stored + store_index].push_back (computed_quantities[0](comp)); + store_index++; + } } } } // end of loop over points @@ -773,9 +792,9 @@ void PointValueHistory 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 ()); - } + { + 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 @@ -783,16 +802,12 @@ void PointValueHistory // 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 + // Repeat for component_mask + typename std::map >::iterator mask = component_mask.find(vector_name); + Assert (mask != component_mask.end(), ExcMessage("vector_name not in class")); + Assert (mask->second.size () == dof_handler->get_fe ().n_components (), ExcDimensionMismatch (mask->second.size (), dof_handler->get_fe ().n_components ())); + unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true); typename std::vector >::iterator point = point_geometry_data.begin (); Vector value (dof_handler->get_fe().n_components()); @@ -801,16 +816,17 @@ void PointValueHistory // 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); + 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++) + // Look up the component_mask and add + // in components according to that mask + for (unsigned int store_index = 0, comp = 0; comp < mask->second.size(); comp++) { - data_store_field->second[data_store_index * components->second.size() + comp].push_back (value (components->second[comp])); + if (mask->second[comp]) + { + data_store_field->second[data_store_index * n_stored + store_index].push_back (value (comp)); + store_index++; + } } } } @@ -975,20 +991,23 @@ void PointValueHistory for (; data_store_begin != data_store.end (); data_store_begin++) { - unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components + typename std::map >::iterator mask = component_mask.find(data_store_begin->first); + unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true); 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] << "> "; - } + AssertThrow (names.size() == n_stored, ExcDimensionMismatch (names.size(), n_stored)); + for (unsigned int component = 0; component < names.size(); component++) + { + to_gnuplot << "<" << names[component] << "> "; + } } else - { - for (unsigned int component = 0; component < n_comps; component++) + { + for (unsigned int component = 0; component < n_stored; component++) { - to_gnuplot << "<" << data_store_begin->first << "_" << component << "> "; + to_gnuplot << "<" << data_store_begin->first << "_" << component << "> "; } } } @@ -1008,11 +1027,12 @@ void PointValueHistory for (; data_store_begin != data_store.end (); data_store_begin++) { - unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components + typename std::map >::iterator mask = component_mask.find(data_store_begin->first); + unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true); - for (unsigned int component = 0; component < n_comps; component++) + for (unsigned int component = 0; component < n_stored; component++) { - to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_comps + component][key]; + to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_stored + component][key]; } } to_gnuplot << "\n"; @@ -1183,25 +1203,29 @@ void PointValueHistory std::map > >::iterator data_store_begin = data_store.begin (); + if (data_store_begin != data_store.end()) + { + out << "Mnemonic: data set size (mask size, n true components) : n data sets\n"; + } 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 mask = component_mask.find(vector_name); + Assert (mask != component_mask.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 () << " ("; - out << components->second.size() << ") : "; + out << mask->second.size() << ", " << std::count (mask->second.begin(), mask->second.end(), true) << ") : "; out << (data_store_begin->second)[0].size () << "\n"; } else { out << data_store_begin->first << ": " << data_store_begin->second.size () << " ("; - out << components->second.size() << ") : "; + out << mask->second.size() << ", " << std::count (mask->second.begin(), mask->second.end(), true) << ") : "; out << "No points added" << "\n"; } // add names, if available @@ -1285,17 +1309,17 @@ template void PointValueHistory ::tria_change_listener () { - // this function is called by the - // Triangulation whenever something - // changes, by virtue of having - // attached the function to the - // signal handler in the - // triangulation object - - // we record the fact that the mesh - // has changed. we need to take - // this into account next time we - // evaluate the solution + // this function is called by the + // Triangulation whenever something + // changes, by virtue of having + // attached the function to the + // signal handler in the + // triangulation object + + // we record the fact that the mesh + // has changed. we need to take + // this into account next time we + // evaluate the solution triangulation_changed = true; } -- 2.39.5