From: rapson Date: Wed, 8 Jul 2009 15:43:50 +0000 (+0000) Subject: Adding point_value_history class to numerics and compute_projection_from_face_quadrat... X-Git-Tag: v8.0.0~7517 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=28a1d738d4ff7e4f9d8b5c2fed1e70e404c19cfd;p=dealii.git Adding point_value_history class to numerics and compute_projection_from_face_quadrature_points_matrix to FETools git-svn-id: https://svn.dealii.org/trunk@19042 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index 9eebc51379..6fae20c4f4 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -644,6 +644,31 @@ class FETools compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, const Quadrature &quadrature, FullMatrix &I_q); + + + + + /** + * This method implements the + * FETools::compute_projection_from_quadrature_points_matrix + * method for faces of a mesh. + * The matrix that it returns, X, is face specific + * and its size is fe.dofs_per_cell by + * rhs_quadrature.size(). + * The dimension, dim must be larger than 1 for this class, + * since Quadrature objects are required. See the + * documentation on the Quadrature class for more information. + */ + template + static + void + compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const typename DoFHandler::active_cell_iterator & cell, + unsigned int face, + FullMatrix &X); + //@} @@ -1330,6 +1355,17 @@ class FETools */ DeclException1(ExcLeastSquaresError, double, << "Least squares fit leaves a gap of " << arg1); + + /** + * Exception thrown if one variable + * may not be greater than another. + * + * @ingroup Exceptions + */ + DeclException2 (ExcNotGreaterThan, + int, int, + << arg1 << " must be greater than " << arg2); + private: /** * Return a finite element that diff --git a/deal.II/deal.II/include/numerics/point_value_history.h b/deal.II/deal.II/include/numerics/point_value_history.h new file mode 100644 index 0000000000..45f392c296 --- /dev/null +++ b/deal.II/deal.II/include/numerics/point_value_history.h @@ -0,0 +1,365 @@ +//--------------------------------------------------------------------------- +// +// Copyright (C) 2009 by Michael Rapson and the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#ifndef __dealii__point_value_history_h +#define __dealii__point_value_history_h + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + namespace PointValueHistory + { + template class PointGeometryData; + } +} + + + +namespace internal +{ + namespace PointValueHistory + { + /*! + * A class that stores the data needed to reference the support points closest to one + * requested point. + */ + template + class PointGeometryData + { + public: + PointGeometryData(std::vector > new_locations, std::vector new_sol_indices); + std::vector > locations; + std::vector solution_indices; + }; + } +} + + + +/*! + * PointValueHistory tackles the overhead of plotting time (or any other iterative process) graphs + * of solution values at specific points on the mesh. The user specifies the points which the solution + * should be monitored at ahead of time, as well as giving each solution vector that they want to + * record a mnemonic name. Then for each + * step the user 'pushes back' the data available from that time step and the class extracts + * data for the requested points to store it. Finally once the computation is finished, the user can + * request Gnuplot output files to be generated. + * + * The user can store extra variables which do not relate to mesh location by the parameter by specifying + * n_independent_variables. The class then expects a std::vector of size n_independent_variables to be added + * during each step using the method @p push_back_independent. This may be used for example for recording + * external input, logging solver performance data such as time taken to solve the step and solver steps + * before convergence, or saving norms calculated. + * + * Currently the code selects the nearest support points to a given point to extract data from. This makes the + * code run at each time step very short, since looping over the mesh to extract the needed dof_index + * can be done just once at the start. However this does lead to problems when FiniteElements which do not + * assign dofs to actual mesh locations are used, (FEs without support points). The FEFieldFunction class + * allows users to evaluate a solution at any point in dof handler and will work even for FEs without support + * points. This functionality is not currently offered through this class since it has more overhead per iteration + * than the approach implemented, but this class can be extended to offer it. + * + * Code snippets showing common usage of the class flows: + * + * @code + * #include + * //.... + * + * //... code to setup Triangulation, perform any refinement necessary + * // and setup DoFHandler, sizing solution Vectors etc + * + * // call the constructor + * unsigned int n_inputs = 1; // just one independent value, which happens to be an input + * PointValueHistory node_monitor(dof_handler, n_inputs); + * // or if node_monitor_2 is defined in class + * node_monitor_2 = PointValueHistory (dof_handler, n_inputs); + * + * // setup fields and points required + * node_monitor.add_field_name("Solution"); + * std::vector > point_vector(2, Point ::Point()); + * point_vector[0] = Point ::Point(0, 0); + * point_vector[1] = Point ::Point(0.25, 0); + * node_monitor.add_points(point_vector); // multiple points at once + * node_monitor.add_point(Point::Point(1, 0.2)); // add a single point + * node_monitor.close(); // close the class once the setup is complete + * node_monitor.status(std::cout); // print out status to check if desired + * + * // ... more code ... + * + * // ... in an iterative loop ... + * // double time, vector with size 1 input_value, + * // and Vector solution calculated in the loop + * node_monitor.start_new_dataset(time); + * node_monitor.push_back_independent(input_value); + * node_monitor.evaluate_field("Solution", solution); + * + * // ... end of iterative loop ... + * + * node_monitor.write_gnuplot("node"); // write out data files + * node_monitor_2.clear (); //clear dof_handler (only if required) + * + * @endcode + */ +template +class PointValueHistory +{ +public: + /*! + * Provide a stripped down instance of the class which does not support adding points or mesh data. + * This may be used for example for recording external input or logging solver performance data. + */ + PointValueHistory (unsigned int n_independent_variables = 0); + + /*! + * Constructor linking the class to a specific @p DoFHandler. Since 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. + * Triangulation 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. + * + * The user can store extra variables which do not relate to mesh location by specifying the + * number required using n_independent_variables and making calls to @p push_back_independent as needed. + * This may be used for example for recording external input or logging solver performance data. + */ + PointValueHistory (const DoFHandler & dof_handler, unsigned int n_independent_variables = 0); + + /*! + * Copy constructor explicitly provided, although default copy operator would be sufficient. + * This constructor can be safely called with a @p PointValueHistory object that contains data, but + * this could be expensive and should be avoided. + */ + PointValueHistory (const PointValueHistory & point_value_history); + + /*! + * Assignment operator explicitly provided, although default assignment operator would be sufficient. + * This assignment operator can be safely called once the class is closed and data added, but this + * is provided primarily to allow a @p PointValueHistory object declared in a class to be reinitialized + * later in the class. Using the assignment operator when the object contains data could be expensive. + */ + PointValueHistory & operator=(const PointValueHistory & point_value_history); + + /*! + Deconstructor, nothing actually required! + */ + ~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 for use when @p evaluate_field is called. If more than one point is required + * rather used the @p add_points method since this minimizes iterations over the mesh. + */ + void add_point(const Point & location); + + /*! + * 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 for use when @p evaluate_field is called. If more than one point is required + * rather call this method as it is more efficient than the add_points method since it minimizes + * iterations over the mesh. The points are added to the internal database in the order they appear in the + * list and there is always a one to one correspondence between the requested point and the added point, + * even if a point is requested multiple times. + */ + void add_points (const std::vector > & locations); + + + + /*! + * Put another mnemonic string (and hence @p VECTOR) into the class. This also adds extra entries for + * points that are already in the class, so @p add_field_name and @p add_points can be called in any order. + */ + void add_field_name(const std::string &vector_name); + + /*! + * Extract values at the stored points from the VECTOR supplied and add them to the new dataset + * in vector_name. If a @p DoFHandler is used, this method must be called for each dataset (time step, + * iteration, etc) for each vector_name, otherwise a @p ExcDataLostSync error can occur. + */ + template + void evaluate_field(const std::string &vector_name, const VECTOR & solution); + + /*! + * Add the key for the current dataset to the dataset. Although calling this method first is + * sensible, the order in which this method, @p evaluate_field and @p push_back_independent is + * not important. It is however important that all the data for a give dataset is added to each + * dataset and that it is added before a new data set is started. This prevents a @p ExcDataLostSync. + */ + void start_new_dataset(double key); + + /*! + * If independent values have been set up, this method stores these values. This should only be called + * once per dataset, and if independent values are used it must be called for every dataset. A @p ExcDataLostSync + * exception can be thrown if this method is not called. + */ + void push_back_independent(const std::vector &independent_values); + + + /*! + * Write out a series of .gpl files named base_name + "-00.gpl", base_name + "-01.gpl" etc. The + * data file gives info about where the support points selected and interpreting the data. + * If @p n_indep != 0 an additional file base_name + "_indep.gpl" containing key and independent + * data. The file name agrees with the order the points were added to the class. + */ + void write_gnuplot(const std::string &base_name); + + + /*! + * Return a @p Vector with the indices of selected points flagged with a 1. This method is mainly + * for testing and verifying that the class is working correctly. By passing this vector to a + * DataOut object, the user can verify that the positions returned by @p PointValueHistory< dim >::get_points + * agree with the positions that @p DataOut interprets from the @p Vector returned. The code snippet below + * demonstrates how this could be done: + * @code + // Make a DataOut object and attach the dof_handler + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + * + // Call the mark_locations method to get the vector with indices flagged + Vector node_locations = node_monitor.mark_locations(); + * + // Add the vector to the data_out object and write out a file in the usual way + data_out.add_data_vector(node_locations, "Monitor_Locations"); + data_out.build_patches(2); + std::ofstream output("locations.gpl"); + data_out.write_gnuplot(output); @endcode + */ + + Vector mark_locations(); + + + /*! + * Stores the actual location of each support point selected by the @p add_point(s) method. + * This can be used to compare with the point requested, for example by using the + * @p Point::distance function. For convenience, location is resized to the + * correct number of points by the method. + */ + void get_points (std::vector > > & locations); + + + /*! + * Once datasets have been added to the class, requests to add additional points + * will make the data interpretation unclear. The boolean @p closed defines a state of the + * class and ensures this does not happen. Additional points or vectors can only be added + * while the class is not closed, and the class must be closed before datasets can be added or + * written to file. @p PointValueHistory::get_points and @p PointValueHistory::status do not require + * the class to be closed. If a method that requires a class to be open or close is called while in + * the wrong state a @p ExcInvalidState exception is thrown. + */ + void close(); + + + /*! + * Delete the lock this object has to the @p DoFHandler used the last time the class was created. + * This method should not normally need to be called, but can be useful to ensure that the + * @p DoFHandler is released before it goes out of scope if the @p PointValueHistory class might + * live longer than it. Once this method has been called, the majority of methods will throw a + * @p ExcInvalidState exception, so if used this method should be the last call to the class. + */ + void clear(); + + /*! + * Print useful debugging information about the class, include details about which support + * points were selected for each point and sizes of the data stored. + */ + void status(std::ostream &out); + + + /*! + * 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 true if all sizes are within 1 of each other (needed to allow data to be added), + * with @p strict = @p true they must be exactly equal. + */ + + bool deep_check (bool strict); + + /*! + * A call has been made to @p push_back_independent when no independent values were requested. + */ + DeclException0(ExcNoIndependent); + + /*! + * This error is thrown to indicate that the data sets appear to be out of sync. + * The class requires that the number of dataset keys is the same as the number of independent + * values sets and mesh linked value sets. The number of each of these is allowed to differ by one to + * allow new values to be added with out restricting the order the user choses to do so. + * Special cases of no @p DoFHandler and no independent values should not trigger this error. + */ + DeclException0(ExcDataLostSync); + + + /*! + * 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 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 that + * the internal dof indices stored are no longer meaningful. + */ + DeclException2(ExcDoFHandlerChanged, + int, + int, + << "Original n_dofs (" << arg1 << ") != current n_dofs (" << arg2 << ")"); + +private: + + std::vector dataset_key; // commonly time, but possibly time step, iteration etc + std::vector > independent_values; // values that do not depend on grid location + std::map > > data_store; // Saves data for each mnemonic entry. data_store: mnemonic -> [component] [time_instance] + std::vector > point_geometry_data; // Saves the location and other mesh info about support points + std::pair > > pair_data; // could possibly be removed, just used as the value to be added to data_store. + bool closed; + bool cleared; + + SmartPointer > 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. + */ + + unsigned int n_dofs; + unsigned int n_indep; + +}; + + +DEAL_II_NAMESPACE_CLOSE +#endif /* __dealii__point_value_history_h */ diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 9e2c36dbd2..dd34c061ca 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -1958,6 +1958,68 @@ compute_interpolation_to_quadrature_points_matrix (const FiniteElement +void +FETools:: +compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const typename DoFHandler::active_cell_iterator & cell, + unsigned int face, + FullMatrix &X) +{ + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree)); + + + + // build the matrices M and Q + // described in the documentation + FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + + { + // need an FEFaceValues object to evaluate shape function + // values on the specified face. + FEFaceValues fe_face_values (fe, lhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + + for (unsigned int i=0; i fe_face_values (fe, rhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); + + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); + M_inverse.mmult (X, Q); + + Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); + Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); +} + + + + /*-------------- Explicit Instantiations -------------------------------*/ @@ -2412,6 +2474,18 @@ compute_interpolation_to_quadrature_points_matrix (const FiniteElement &quadrature, FullMatrix &I_q); +#if deal_II_dimension != 1 +template +void +FETools:: +compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const DoFHandler::active_cell_iterator & cell, + unsigned int face, + FullMatrix &X); +#endif + /*---------------------------- fe_tools.cc ---------------------------*/ diff --git a/deal.II/deal.II/source/numerics/point_value_history.cc b/deal.II/deal.II/source/numerics/point_value_history.cc new file mode 100644 index 0000000000..129bcd5ec7 --- /dev/null +++ b/deal.II/deal.II/source/numerics/point_value_history.cc @@ -0,0 +1,708 @@ +//--------------------------------------------------------------------------- +// +// Copyright (C) 2009 by Michael Rapson and the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include + + + +DEAL_II_NAMESPACE_OPEN + +#include "point_value_history.inst" + +/// Only a constructor needed for this class (a struct really) +template +internal::PointValueHistory::PointGeometryData +::PointGeometryData (std::vector > new_locations, std::vector new_sol_indices) +{ + locations = new_locations; + solution_indices = new_sol_indices; +} + + + +template +PointValueHistory +::PointValueHistory (unsigned int n_independent_variables) : +dof_handler (0, "No DoFHandler"), +n_dofs (0), +n_indep (n_independent_variables) +{ + closed = false; + cleared = false; + // make a vector for keys + dataset_key = std::vector (); // initialize the std::vector + + // make a vector of independent values + independent_values = std::vector > (n_indep, std::vector (0)); +} + + + +template +PointValueHistory +::PointValueHistory (const DoFHandler & dof_handler, unsigned int n_independent_variables) : +dof_handler (&dof_handler, typeid (*this).name ()), +n_dofs (dof_handler.n_dofs ()), +n_indep (n_independent_variables) +{ + closed = false; + cleared = false; + // make a vector to store keys + dataset_key = std::vector (); // initialize the std::vector + + // make a vector for the independent values + independent_values = std::vector > (n_indep, std::vector (0)); +} + + + +template +PointValueHistory +::PointValueHistory (const PointValueHistory & point_value_history) +{ + dataset_key = point_value_history.dataset_key; + independent_values = point_value_history.independent_values; + data_store = point_value_history.data_store; + 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; + n_indep = point_value_history.n_indep; +} + + + +template + PointValueHistory & PointValueHistory::operator= (const PointValueHistory & point_value_history) +{ + dataset_key = point_value_history.dataset_key; + independent_values = point_value_history.independent_values; + data_store = point_value_history.data_store; + 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; + n_indep = point_value_history.n_indep; + return * this; +} + + + +template +PointValueHistory +::~PointValueHistory () +{ +} + + + +template +void PointValueHistory +::add_point (const Point & location) +{ + // 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 (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ()); + + // 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. + 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 + 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 + 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); + std::vector > current_points (n_components, Point ()); + for (unsigned int support_point = 0; support_point < n_support_points; support_point++) + { + // setup valid data in the empty vectors + unsigned int component = 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 + for (; cell != endc; cell++) + { + fe_values.reinit (cell); + + for (unsigned int support_point = 0; support_point < n_support_points; support_point++) + { + unsigned int component = dof_handler->get_fe ().system_to_component_index (support_point).first; + Point test_point = fe_values.quadrature_point (support_point); + + if (location.distance (test_point) < location.distance (current_points [component])) + { + // save the data + current_points [component] = test_point; + current_cell = cell; + current_fe_index [component] = support_point; + } + } + } + + + std::vector 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. + + + 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]]); + } + + internal::PointValueHistory::PointGeometryData new_point_geometry_data (current_points, new_solution_indices); + point_geometry_data.push_back (new_point_geometry_data); + + std::map > >::iterator + data_store_begin = data_store.begin (); + for (; data_store_begin != data_store.end (); data_store_begin++) + { + // add an extra row to each vector entry + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + data_store_begin->second.push_back (std::vector (0)); + } + } +} + + + +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! + + // 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 (dof_handler->get_fe ().has_support_points (), ExcNotImplemented ()); + + // 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. + 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 + 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 + + // need to store one value per point per component + std::vector ::active_cell_iterator > current_cell (locations.size (), cell); + + fe_values.reinit (cell); + std::vector > temp_points (n_components, Point ()); + 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 + unsigned int component = dof_handler->get_fe ().system_to_component_index (support_point).first; + temp_points [component] = fe_values.quadrature_point (support_point); + temp_fe_index [component] = support_point; + } + std::vector > > 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 + for (; cell != endc; cell++) + { + fe_values.reinit (cell); + for (unsigned int support_point = 0; support_point < n_support_points; support_point++) + { + unsigned int component = dof_handler->get_fe ().system_to_component_index (support_point).first; + Point test_point = fe_values.quadrature_point (support_point); + + for (unsigned int point = 0; point < locations.size (); point++) + { + if (locations[point].distance (test_point) < locations[point].distance (current_points[point][component])) + { + // save the data + current_points[point][component] = test_point; + current_cell[point] = cell; + current_fe_index[point][component] = support_point; + } + } + } + } + + std::vector local_dof_indices (dof_handler->get_fe ().dofs_per_cell); + for (unsigned int point = 0; point < locations.size (); point++) + { + current_cell[point]->get_dof_indices (local_dof_indices); + std::vector new_solution_indices; + + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + new_solution_indices.push_back (local_dof_indices[current_fe_index[point][component]]); + } + + internal::PointValueHistory::PointGeometryData new_point_geometry_data (current_points[point], new_solution_indices); + + point_geometry_data.push_back (new_point_geometry_data); + + std::map > >::iterator + data_store_begin = data_store.begin (); + for (; data_store_begin != data_store.end (); data_store_begin++) + { + // add an extra row to each vector entry + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + data_store_begin->second.push_back (std::vector (0)); + } + } + } +} + + + + + + +template +void PointValueHistory +::add_field_name (const std::string &vector_name) +{ + // can't be closed to add additional points or vectors + AssertThrow (!closed, ExcInvalidState ()); + AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); + AssertThrow (!cleared, ExcInvalidState ()); + AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + + + // make and add a new vector point_geometry_data.size() long + pair_data.first = vector_name; + int n_datastreams = point_geometry_data.size () * (dof_handler->get_fe ().n_components ()); // each point has n_components sub parts + std::vector < std::vector > vector_size (n_datastreams, + std::vector (0)); + pair_data.second = vector_size; + data_store.insert (pair_data); +} + + + +template +void PointValueHistory +::close () +{ + closed = true; +} + + + +template +void PointValueHistory +::clear () +{ + cleared = true; + dof_handler = SmartPointer > (0, "Cleared"); +} + +// Need to test that the internal data has a full and complete dataset for each key. That is that +// the data has not got 'out of sync'. Testing that dataset_key is within 1 of independent_values +// is cheap and is done in all three methods. Evaluate_field will check that its vector_name is within +// 1 of dataset_key. However this leaves the possibility that the user has neglected to call evaluate_field +// on one vector_name consistently. To catch this case start_new_dataset will call bool deap_check () which +// will test all vector_names and return a bool. This can be called from an Assert statement. + + + +template +template +void PointValueHistory +::evaluate_field (const std::string &vector_name, const VECTOR & solution) +{ + // 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 ())); + 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 ()); + } + typename std::vector >::iterator node = point_geometry_data.begin (); + for (unsigned int data_store_index = 0; node != point_geometry_data.end (); node++, data_store_index++) + { + // step through each node, to access the Solution_index and the vector index + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + unsigned int solution_index = node->solution_indices[component]; + (data_store[vector_name])[data_store_index * dof_handler->get_fe ().n_components () + component].push_back (solution (solution_index)); + } + } +} + + + +template +void PointValueHistory +::start_new_dataset (double key) +{ + // must be closed to add data to internal members. + Assert (closed, ExcInvalidState ()); + Assert (!cleared, ExcInvalidState ()); + Assert (deep_check (false), ExcDataLostSync ()); + + dataset_key.push_back (key); +} + + + +template +void PointValueHistory +::push_back_independent (const std::vector &indep_values) +{ + // 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)); + Assert (n_indep != 0, ExcNoIndependent ()); + Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ()); + + for (unsigned int component = 0; component < n_indep; component++) + independent_values[component].push_back (indep_values[component]); +} + + + +template +void PointValueHistory +::write_gnuplot (const std::string &base_name) +{ + AssertThrow (closed, ExcInvalidState ()); + AssertThrow (!cleared, ExcInvalidState ()); + AssertThrow (deep_check (true), ExcDataLostSync ()); + + // write inputs to a file + if (n_indep != 0) + { + std::string filename = base_name + "_indep.gpl"; + std::ofstream to_gnuplot (filename.c_str ()); + + to_gnuplot << "# Data independent of mesh location\n"; + + // write column headings + std::map > >::iterator + data_store_begin = data_store.begin (); + to_gnuplot << "# "; + + 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++) + { + std::map > >::iterator + data_store_begin = data_store.begin (); + to_gnuplot << dataset_key[key]; + + for (unsigned int component = 0; component < n_indep; component++) + { + to_gnuplot << " " << independent_values[component][key]; + } + to_gnuplot << "\n"; + } + + to_gnuplot.close (); + } + + + + // write points to a file + if (n_dofs != 0) + { + AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); + AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + + unsigned int n_components = dof_handler->get_fe ().n_components (); + + typename std::vector >::iterator node = point_geometry_data.begin (); + for (unsigned int data_store_index = 0; node != point_geometry_data.end (); node++, data_store_index++) + { + // 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 upto 100 points + std::ofstream to_gnuplot (filename.c_str ()); + + // put helpful info about the support point into the file as comments + to_gnuplot << "# DoF_index : Location (for each component)\n"; + for (unsigned int component = 0; component < n_components; component++) + { + to_gnuplot << "# " << node->solution_indices[component] << " : " << node->locations [component] << "\n"; + } + to_gnuplot << "\n"; + + // write column headings + std::map > >::iterator + data_store_begin = data_store.begin (); + to_gnuplot << "# "; + + for (unsigned int component = 0; component < n_indep; component++) + { + to_gnuplot << " "; + } + + for (; data_store_begin != data_store.end (); data_store_begin++) + { + for (unsigned int component = 0; component < n_components; component++) + { + to_gnuplot << "<" << data_store_begin->first << "_" << component << "> "; + } + } + to_gnuplot << "\n"; + + // write data stored for the node + for (unsigned int key = 0; key < dataset_key.size (); key++) + { + std::map > >::iterator + data_store_begin = data_store.begin (); + to_gnuplot << dataset_key[key]; + + for (unsigned int component = 0; component < n_indep; component++) + { + to_gnuplot << " " << independent_values[component][key]; + } + + for (; data_store_begin != data_store.end (); data_store_begin++) + { + for (unsigned int component = 0; component < n_components; component++) + { + to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_components + component][key]; + } + } + to_gnuplot << "\n"; + } + + to_gnuplot.close (); + } + } +} + + + +template +Vector PointValueHistory +::mark_locations () +{ + // a method to put a one at each point on the grid where a location is defined + AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); + AssertThrow (!cleared, ExcInvalidState ()); + AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + + Vector dof_vector (dof_handler->n_dofs ()); + + typename std::vector >::iterator node = point_geometry_data.begin (); + for (; node != point_geometry_data.end (); node++) + { + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + dof_vector (node->solution_indices[component]) = 1; + } + } + return dof_vector; +} + + + +template +void PointValueHistory +::get_points (std::vector > > & locations) +{ + AssertThrow (n_dofs != 0, ExcDoFHandlerRequired ()); + AssertThrow (!cleared, ExcInvalidState ()); + AssertThrow (n_dofs == dof_handler->n_dofs (), ExcDoFHandlerChanged (n_dofs, dof_handler->n_dofs ())); + + std::vector > > actual_points; + typename std::vector >::iterator node = point_geometry_data.begin (); + + for (; node != point_geometry_data.end (); node++) + { + actual_points.push_back (node->locations); + } + locations = actual_points; +} + + + +template +void PointValueHistory +::status (std::ostream &out) +{ + out << "***PointValueHistory status output***\n\n"; + out << "Closed: " << closed << "\n"; + out << "Cleared: " << cleared << "\n"; + out << "Geometric Data" << "\n"; + + typename std::vector >::iterator node = point_geometry_data.begin (); + if (node == point_geometry_data.end ()) + { + out << "No nodes stored currently\n"; + } + else + { + if (!cleared) + { + out << "# DoF_index : Location (for each component)\n"; + for (; node != point_geometry_data.end (); node++) + { + for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++) + { + out << node->solution_indices[component] << " : " << node->locations [component] << "\n"; + } + out << "\n"; + } + } + else + { + out << "#Cannot access DoF_indices once cleared\n"; + } + } + out << "\n"; + + std::cout << dataset_key.size () << " data keys are stored \n"; + + + if (independent_values.size () != 0) + { + out << "Independent value(s): " << independent_values.size () << " : " << independent_values[0].size () << "\n"; + } + else + { + out << "No independent values stored\n"; + } + + std::map > >::iterator + data_store_begin = data_store.begin (); + for (; data_store_begin != data_store.end (); data_store_begin++) + { + if (data_store_begin->second.size () != 0) + { + out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << (data_store_begin->second)[0].size () << "\n"; + } + else + { + out << data_store_begin->first << ": " << data_store_begin->second.size () << " : " << "No points added" << "\n"; + } + } + out << "\n"; + out << "***end of status output***\n\n"; +} + + + +template +bool PointValueHistory +::deep_check (bool strict) +{ + // test ways that it can fail, if control reaches last statement return true + if (strict) + { + if (n_indep != 0) + { + if (dataset_key.size () != independent_values[0].size ()) + { + return false; + } + } + std::map > >::iterator + data_store_begin = data_store.begin (); + if (n_dofs != 0) + { + for (; data_store_begin != data_store.end (); data_store_begin++) + { + if ((data_store_begin->second)[0].size () != dataset_key.size ()) + return false; + // this loop only tests one member for each name, i.e. checks the user + // it will not catch internal errors which do not update all fields for a name. + } + } + return true; + } + if (n_indep != 0) + { + if (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) >= 2) + { + return false; + } + } + + if (n_dofs != 0) + { + std::map > >::iterator + 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. + } + } + return true; +} + + +// explicit instantiations +template class PointValueHistory; + + +DEAL_II_NAMESPACE_CLOSE + diff --git a/deal.II/deal.II/source/numerics/point_value_history.inst.in b/deal.II/deal.II/source/numerics/point_value_history.inst.in new file mode 100644 index 0000000000..77fd01c549 --- /dev/null +++ b/deal.II/deal.II/source/numerics/point_value_history.inst.in @@ -0,0 +1,7 @@ +for (VEC : SERIAL_VECTORS) +{ + template + void PointValueHistory::evaluate_field + (const std::string &, + const VEC &); +}