From 4a92c9c4c73f94673f5fde546f5dbaf7f9d2b354 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 19 Jul 1999 07:57:48 +0000 Subject: [PATCH] Take these files over from the wave project. git-svn-id: https://svn.dealii.org/trunk@1594 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/numerics/data_out_stack.h | 400 ++++++++++++++++++ .../deal.II/source/numerics/data_out_stack.cc | 332 +++++++++++++++ 2 files changed, 732 insertions(+) create mode 100644 deal.II/deal.II/include/numerics/data_out_stack.h create mode 100644 deal.II/deal.II/source/numerics/data_out_stack.cc diff --git a/deal.II/deal.II/include/numerics/data_out_stack.h b/deal.II/deal.II/include/numerics/data_out_stack.h new file mode 100644 index 0000000000..b671b75376 --- /dev/null +++ b/deal.II/deal.II/include/numerics/data_out_stack.h @@ -0,0 +1,400 @@ +/*---------------------------- data_out_stack.h ---------------------------*/ +/* $Id$ */ +#ifndef __data_out_stack_H +#define __data_out_stack_H +/*---------------------------- data_out_stack.h ---------------------------*/ + + +#include +#include +#include +#include + +#include +#include + + + +/** + * This class is used to stack the output from several computations + * into one output file by stacking the data sets in another + * co-ordinate direction orthogonal to the space directions. The most + * common use is to stack the results of several time steps into one + * space-time output file, or for example to connect the results of + * solutions of a parameter dependent equation for several parameter + * value together into one. The interface is mostly modelled after the + * #DataOut# class. + * + * We will explain the concept for a time dependent problem, but + * instead of the time any parameter can be substituted. In our + * example, a solution of an equation is computed for each discrete + * time level. This is then added to an object of the present class + * and after all time levels are added, a space-time plot will be + * written in any of the output formats supported by the base + * class. Upon output, the (spatial) solution on each time level is + * extended into the time direction by writing it twice, once for the + * time level itself and once for a time equal to the time level minus + * a given time step. These two copies are connected, to form a + * space-time slab, with constant values in time. + * + * Due to the piecewise constant output in time, the written solution + * will in general be discontinuous at discrete time levels, but the + * output is still sufficient in most cases. More sophisticated + * interpolations in time may be added in the future. + * + * + * \secion{Example of Use} + * + * The following little example shall illustrate the different steps + * of use of this class. It is assumed that the finite element used is + * composed of two components, #u# and #v#, that the solution vector + * is named #solution# and that a vector #error# is computed which + * contains an error indicator for each spatial cell. + * + * Note that unlike for the #DataOut# class it is necessary to first + * declare data vectors and the names of the components before first + * use. This is because on all time levels the same data should be + * present to produce reasonable time-space output. The output is + * generated with two subdivisions in each space and time direction, + * which is suitable for quadratic finite elements in space, for + * example. + * + * \begin{verbatim} + * DataOutStack data_out_stack; + * + * // first declare the vectors + * // to be used later + * vector solution_names; + * solution_names.push_back ("u"); + * solution_names.push_back ("v"); + * data_out_stack.declare_data_vector (solution_names, + * DataOutStack::dof_data); + * data_out_stack.declare_data_vector ("error", + * DataOutStack::cell_data); + * + * // now do computations + * for (double parameter=0; ...) + * { + * DoFHandler dof_handler; + * ... // compute something + * + * // now for output + * data_out_stack.new_parameter_value (parameter, + * delta_parameter); + * data_out_stack.attach_dof_handler (dof_handler); + * data_out_stack.add_data_vector (solution, solution_names); + * data_out_stack.add_data_vector (error, "error"); + * data_out_stack.build_patches (2); + * data_out_stack.finish_parameter_value (); + * }; + * \end{verbatim} + * + * @author Wolfgang Bangerth, 1999 + */ +template +class DataOutStack : public DataOutInterface +{ + public: + /** + * Data type declaring the two types + * of vectors which are used in this + * class. + */ + enum VectorType { cell_vector, dof_vector }; + + /** + * Destructor. Only declared to make + * it #virtual#. + */ + virtual ~DataOutStack (); + + /** + * Start the next set of data for a + * specific parameter value. The argument + * #parameter_step# denotes the interval + * (in backward direction, counted from + * #parameter_value#) with which the + * output will be extended in parameter + * direction, i.e. orthogonal to the + * space directions. + */ + void new_parameter_value (const double parameter_value, + const double parameter_step); + + /** + * Attach the DoF handler for the + * grid and data associated with the + * parameter previously set by + * #new_parameter_value#. + * + * This has to happen before adding + * data vectors for the present + * parameter value. + */ + void attach_dof_handler (const DoFHandler &dof_handler); + + /** + * Declare a data vector. The #vector_type# + * argument determines whether the data + * vector will be considered as DoF or + * cell data. + * + * This version may be called if the + * finite element presently used by the + * #DoFHandler# (and previously attached + * to this object) has only one component + * and therefore only one name needs to + * be given. + */ + void declare_data_vector (const string &name, + const VectorType vector_type); + + /** + * Declare a data vector. The #vector_type# + * argument determines whether the data + * vector will be considered as DoF or + * cell data. + * + * This version must be called if the + * finite element presently used by the + * #DoFHandler# (and previously attached + * to this object) has more than one + * component and therefore more than one + * name needs to be given. However, you + * can also call this function with a + * #vector# containing only one + * element if the finite element has + * only one component. + */ + void declare_data_vector (const vector &name, + const VectorType vector_type); + + + /** + * Add a data vector for the presently + * set value of the parameter. + * + * This version may be called if the + * finite element presently used by the + * #DoFHandler# (and previously attached + * to this object) has only one component + * and therefore only one name needs to + * be given. + * + * The data vector must have been + * registered using the #declare_data_vector# + * function before actually using it the + * first time. + * + * Note that a copy of this vector is + * stored until #finish_parameter_value# + * is called the next time, so if you are + * short of memory you may want to call + * this function only after all + * computations involving large matrices + * are already done. + */ + template + void add_data_vector (const Vector &vec, + const string &name); + + /** + * Add a data vector for the presently + * set value of the parameter. + * + * This version must be called if the + * finite element presently used by the + * #DoFHandler# (and previously attached + * to this object) has more than one + * component and therefore more than one + * name needs to be given. However, you + * can also call this function with a + * #vector# containing only one + * element if the finite element has + * only one component. + * + * The data vector must have been + * registered using the #declare_data_vector# + * function before actually using it the + * first time. + * + * Note that a copy of this vector is + * stored until #finish_parameter_value# + * is called the next time, so if you are + * short of memory you may want to call + * this function only after all + * computations involving large matrices + * are already done. + */ + template + void add_data_vector (const Vector &vec, + const vector &names); + + /** + * Actually build the patches for output + * by the base classes from the data + * stored in the data vectors and using + * the previously attached #DoFHandler# + * object. + * + * By #n_subdivisions# you can decide + * into how many subdivisions (in each + * space and parameter direction) each + * patch is divided. This is useful + * if higher order elements are used. + * Note however, that the number of + * subdivisions in parameter direction + * is always the same as the one is space + * direction for technical reasons. + */ + void build_patches (const unsigned int n_subdivisions = 1); + + /** + * Release all data that is no + * more needed once #build_patches# + * was called and all other transactions + * for a given parameter value are done. + * + * Couterpart of #new_parameter_value#. + */ + void finish_parameter_value (); + + /** + * Clear all data presently stored + * in this object. + */ + void clear (); + + /** + * Exception + */ + DeclException0 (ExcNoDoFHandlerSelected); + /** + * Exception + */ + DeclException3 (ExcInvalidVectorSize, + int, int, int, + << "The vector has size " << arg1 + << " but the DoFHandler objects says there are " << arg2 + << " degrees of freedom and there are " << arg3 + << " active cells."); + /** + * Exception + */ + DeclException1 (ExcInvalidCharacter, + string, + << "Please use only the characters [a-zA-Z0-9_<>()] for" << endl + << "description strings since AVS will only accept these." << endl + << "The string you gave was <" << arg1 << ">."); + /** + * Exception + */ + DeclException2 (ExcInvalidNumberOfNames, + int, int, + << "You have to give one name per component in your " + << "data vector. The number you gave was " << arg1 + << ", but the number of components is " << arg2); + /** + * Exception + */ + DeclException1 (ExcVectorNotDeclared, + string, + << "The data vector for which the first component has the name " + << arg1 << " has not been declared before."); + /** + * Exception + */ + DeclException0 (ExcDataNotCleared); + /** + * Exception + */ + DeclException0 (ExcDataAlreadyAdded); + /** + * Exception + */ + DeclException1 (ExcNameAlreadyUsed, + string, + << "You tried to declare a component of a data vector with " + << "the name <" << arg1 << ">, but that name is already used."); + + private: + /** + * Present parameter value. + */ + double parameter; + + /** + * Present parameter step, i.e. + * length of the parameter interval + * to be written next. + */ + double parameter_step; + + /** + * DoF handler to be used for the data + * corresponding to the present parameter + * value. + */ + SmartPointer > dof_handler; + + /** + * List of patches of all past and + * present parameter value data sets. + */ + vector > patches; + + /** + * Structure holding data vectors + * (cell and dof data) for the present + * parameter value. + */ + struct DataVector + { + /** + * Data vector. + */ + Vector data; + + /** + * Names of the different components + * within each such data set. + */ + vector names; + }; + + /** + * List of DoF data vectors. + */ + vector dof_data; + + /** + * List of cell data vectors. + */ + vector cell_data; + + /** + * This is the function through + * which derived classes propagate + * preprocessed data in the form of + * #Patch# structures (declared in + * the base class #DataOutBase#) to + * the actual output function. + */ + virtual const vector > & get_patches () const; + + + /** + * Virtual function through + * which the names of data sets are + * obtained by the output functions + * of the base class. + */ + virtual vector get_dataset_names () const; +}; + + + +/*---------------------------- data_out_stack.h ---------------------------*/ +/* end of #ifndef __data_out_stack_H */ +#endif +/*---------------------------- data_out_stack.h ---------------------------*/ diff --git a/deal.II/deal.II/source/numerics/data_out_stack.cc b/deal.II/deal.II/source/numerics/data_out_stack.cc new file mode 100644 index 0000000000..54307a4531 --- /dev/null +++ b/deal.II/deal.II/source/numerics/data_out_stack.cc @@ -0,0 +1,332 @@ +/* $Id$ */ + +#include "../include/data_out_stack.h" +#include +#include +#include +#include +#include +#include +#include + + +template +DataOutStack::~DataOutStack () +{}; + + + +template +void DataOutStack::new_parameter_value (const double p, + const double dp) +{ + parameter = p; + parameter_step = dp; + + // check whether the user called #finish_...# + // at the end of the previous parameter step + // + // this is to prevent serious waste of + // memory + for (vector::const_iterator i=dof_data.begin(); i!=dof_data.end(); ++i) + Assert (i->data.size() == 0, + ExcDataNotCleared ()); + for (vector::const_iterator i=cell_data.begin(); i!=cell_data.end(); ++i) + Assert (i->data.size() == 0, + ExcDataNotCleared ()); + +}; + + +template +void DataOutStack::attach_dof_handler (const DoFHandler &dof) +{ + dof_handler = &dof; +}; + + + +template +void DataOutStack::declare_data_vector (const string &name, + const VectorType vector_type) +{ + vector names; + names.push_back (name); + declare_data_vector (names, vector_type); +}; + + + +template +void DataOutStack::declare_data_vector (const vector &names, + const VectorType vector_type) +{ + // make sure this function is + // not called after some parameter + // values have already been + // processed + Assert (patches.size() == 0, ExcDataAlreadyAdded()); + + // also make sure that no name is + // used twice + for (vector::const_iterator name=names.begin(); name!=names.end(); ++name) + { + for (vector::const_iterator data_set=dof_data.begin(); + data_set!=dof_data.end(); ++data_set) + for (unsigned int i=0; inames.size(); ++i) + Assert (*name != data_set->names[i], ExcNameAlreadyUsed(*name)); + + for (vector::const_iterator data_set=cell_data.begin(); + data_set!=cell_data.end(); ++data_set) + for (unsigned int i=0; inames.size(); ++i) + Assert (*name != data_set->names[i], ExcNameAlreadyUsed(*name)); + }; + + switch (vector_type) + { + case dof_vector: + dof_data.push_back (DataVector()); + dof_data.back().names = names; + break; + + case cell_vector: + cell_data.push_back (DataVector()); + cell_data.back().names = names; + break; + }; +}; + + +template +template +void DataOutStack::add_data_vector (const Vector &vec, + const string &name) +{ + vector names; + names.push_back (name); + add_data_vector (vec, names); +}; + + + + +template +template +void DataOutStack::add_data_vector (const Vector &vec, + const vector &names) +{ + Assert (dof_handler != 0, ExcNoDoFHandlerSelected ()); + // either cell data and one name, + // or dof data and n_components names + Assert (((vec.size() == dof_handler->get_tria().n_active_cells()) && + (names.size() == 1)) + || + ((vec.size() == dof_handler->n_dofs()) && + (names.size() == dof_handler->get_fe().n_components)), + ExcInvalidNumberOfNames (names.size(), dof_handler->get_fe().n_components)); + for (unsigned int i=0; i()") == string::npos, + ExcInvalidCharacter (names[i])); + + if (vec.size() == dof_handler->n_dofs()) + { + vector::iterator data_vector=dof_data.begin(); + for (; data_vector!=dof_data.end(); ++data_vector) + if (data_vector->names == names) + { + data_vector->data.reinit (vec.size()); + copy (vec.begin(), vec.end(), + data_vector->data.begin()); + break; + }; + Assert (data_vector != dof_data.end(), + ExcVectorNotDeclared (names[0])); + } + else + { + vector::iterator data_vector=cell_data.begin(); + for (; data_vector!=cell_data.end(); ++data_vector) + if (data_vector->names == names) + { + data_vector->data.reinit (vec.size()); + copy (vec.begin(), vec.end(), + data_vector->data.begin()); + break; + }; + Assert (data_vector != cell_data.end(), + ExcVectorNotDeclared (names[0])); + }; +}; + + + +template +void DataOutStack::build_patches (const unsigned int n_subdivisions) +{ + // this is mostly copied from the + // DataOut class + + Assert (dof_handler != 0, ExcNoDoFHandlerSelected()); + + const unsigned int n_components = dof_handler->get_fe().n_components; + const unsigned int n_datasets = dof_data.size() * n_components + + cell_data.size(); + + // first count the cells we want to + // create patches of and make sure + // there is enough memory for that + unsigned int n_patches = 0; + for (DoFHandler::active_cell_iterator cell=dof_handler->begin_active(); + cell != dof_handler->end(); ++cell) + ++n_patches; + + + // before we start the loop: + // create a quadrature rule that + // actually has the points on this + // patch, and an object that + // extracts the data on each + // cell to these points + QTrapez<1> q_trapez; + QIterated patch_points (q_trapez, n_subdivisions); + FEValues fe_patch_values (dof_handler->get_fe(), + patch_points, + update_default); + const unsigned int n_q_points = patch_points.n_quadrature_points; + vector patch_values (n_q_points); + vector > patch_values_system (n_q_points, + Vector(n_components)); + + // add the required number of patches + DataOutBase::Patch default_patch; + default_patch.n_subdivisions = n_subdivisions; + default_patch.data.reinit (n_datasets, n_q_points*(n_subdivisions+1)); + patches.insert (patches.end(), n_patches, default_patch); + + // now loop over all cells and + // actually create the patches + vector >::iterator patch = &patches[patches.size()-n_patches]; + unsigned int cell_number = 0; + for (DoFHandler::active_cell_iterator cell=dof_handler->begin_active(); + cell != dof_handler->end(); ++cell, ++patch, ++cell_number) + { + Assert (patch != patches.end(), ExcInternalError()); + + // first fill in the vertices of the patch + switch (dim) + { + case 1: + patch->vertices[0] = Point(cell->vertex(0)(0), + parameter-parameter_step); + patch->vertices[3] = Point(cell->vertex(0)(0), + parameter); + patch->vertices[1] = Point(cell->vertex(1)(0), + parameter-parameter_step); + patch->vertices[2] = Point(cell->vertex(1)(0), + parameter); + break; + + default: + Assert (false, ExcNotImplemented()); + }; + + + // now fill in the the data values. + // note that the required order is + // with highest coordinate running + // fastest, we need to enter each + // value (n_subdivisions+1) times + // in succession + if (n_datasets > 0) + { + fe_patch_values.reinit (cell); + + // first fill dof_data + for (unsigned int dataset=0; datasetdata(dataset,q*(n_subdivisions+1)+i) = patch_values[q]; + } + else + // system of components + { + fe_patch_values.get_function_values (dof_data[dataset].data, + patch_values_system); + for (unsigned int component=0; componentdata(dataset*n_components+component, + q*(n_subdivisions+1)+i) + = patch_values_system[q](component); + }; + }; + + // then do the cell data + for (unsigned int dataset=0; datasetdata(dataset+dof_data.size()*n_components, + q*(n_subdivisions+1)+i) = value; + }; + }; + }; +}; + + + +template +void DataOutStack::finish_parameter_value () +{ + // release lock on dof handler + dof_handler = 0; + for (vector::iterator i=dof_data.begin(); i!=dof_data.end(); ++i) + i->data.reinit (0); + + for (vector::iterator i=cell_data.begin(); i!=cell_data.end(); ++i) + i->data.reinit (0); +}; + + + +template +const vector > & +DataOutStack::get_patches () const +{ + return patches; +}; + + + +template +vector DataOutStack::get_dataset_names () const +{ + vector names; + for (vector::const_iterator dataset=dof_data.begin(); + dataset!=dof_data.end(); ++dataset) + names.insert (names.end(), dataset->names.begin(), dataset->names.end()); + for (vector::const_iterator dataset=cell_data.begin(); + dataset!=cell_data.end(); ++dataset) + names.insert (names.end(), dataset->names.begin(), dataset->names.end()); + + return names; +}; + + + + +// explicit instantiations +template class DataOutStack; +template void DataOutStack::add_data_vector (const Vector &vec, + const string &name); +template void DataOutStack::add_data_vector (const Vector &vec, + const string &name); + -- 2.39.5