--- /dev/null
+/*---------------------------- data_out_stack.h ---------------------------*/
+/* $Id$ */
+#ifndef __data_out_stack_H
+#define __data_out_stack_H
+/*---------------------------- data_out_stack.h ---------------------------*/
+
+
+#include <base/data_out_base.h>
+#include <lac/forward-declarations.h>
+#include <basic/forward-declarations.h>
+#include <base/smartpointer.h>
+
+#include <string>
+#include <vector>
+
+
+
+/**
+ * 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<dim> data_out_stack;
+ *
+ * // first declare the vectors
+ * // to be used later
+ * vector<string> solution_names;
+ * solution_names.push_back ("u");
+ * solution_names.push_back ("v");
+ * data_out_stack.declare_data_vector (solution_names,
+ * DataOutStack<dim>::dof_data);
+ * data_out_stack.declare_data_vector ("error",
+ * DataOutStack<dim>::cell_data);
+ *
+ * // now do computations
+ * for (double parameter=0; ...)
+ * {
+ * DoFHandler<dim> 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 <int dim>
+class DataOutStack : public DataOutInterface<dim+1>
+{
+ 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<dim> &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<string># containing only one
+ * element if the finite element has
+ * only one component.
+ */
+ void declare_data_vector (const vector<string> &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 <typename number>
+ void add_data_vector (const Vector<number> &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<string># 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 <typename number>
+ void add_data_vector (const Vector<number> &vec,
+ const vector<string> &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<const DoFHandler<dim> > dof_handler;
+
+ /**
+ * List of patches of all past and
+ * present parameter value data sets.
+ */
+ vector<DataOutBase::Patch<dim+1> > patches;
+
+ /**
+ * Structure holding data vectors
+ * (cell and dof data) for the present
+ * parameter value.
+ */
+ struct DataVector
+ {
+ /**
+ * Data vector.
+ */
+ Vector<double> data;
+
+ /**
+ * Names of the different components
+ * within each such data set.
+ */
+ vector<string> names;
+ };
+
+ /**
+ * List of DoF data vectors.
+ */
+ vector<DataVector> dof_data;
+
+ /**
+ * List of cell data vectors.
+ */
+ vector<DataVector> 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<DataOutBase::Patch<dim+1> > & get_patches () const;
+
+
+ /**
+ * Virtual function through
+ * which the names of data sets are
+ * obtained by the output functions
+ * of the base class.
+ */
+ virtual vector<string> get_dataset_names () const;
+};
+
+
+
+/*---------------------------- data_out_stack.h ---------------------------*/
+/* end of #ifndef __data_out_stack_H */
+#endif
+/*---------------------------- data_out_stack.h ---------------------------*/
--- /dev/null
+/* $Id$ */
+
+#include "../include/data_out_stack.h"
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/dof.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <fe/fe.h>
+#include <fe/fe_values.h>
+
+
+template <int dim>
+DataOutStack<dim>::~DataOutStack ()
+{};
+
+
+
+template <int dim>
+void DataOutStack<dim>::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<DataVector>::const_iterator i=dof_data.begin(); i!=dof_data.end(); ++i)
+ Assert (i->data.size() == 0,
+ ExcDataNotCleared ());
+ for (vector<DataVector>::const_iterator i=cell_data.begin(); i!=cell_data.end(); ++i)
+ Assert (i->data.size() == 0,
+ ExcDataNotCleared ());
+
+};
+
+
+template <int dim>
+void DataOutStack<dim>::attach_dof_handler (const DoFHandler<dim> &dof)
+{
+ dof_handler = &dof;
+};
+
+
+
+template <int dim>
+void DataOutStack<dim>::declare_data_vector (const string &name,
+ const VectorType vector_type)
+{
+ vector<string> names;
+ names.push_back (name);
+ declare_data_vector (names, vector_type);
+};
+
+
+
+template <int dim>
+void DataOutStack<dim>::declare_data_vector (const vector<string> &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<string>::const_iterator name=names.begin(); name!=names.end(); ++name)
+ {
+ for (vector<DataVector>::const_iterator data_set=dof_data.begin();
+ data_set!=dof_data.end(); ++data_set)
+ for (unsigned int i=0; i<data_set->names.size(); ++i)
+ Assert (*name != data_set->names[i], ExcNameAlreadyUsed(*name));
+
+ for (vector<DataVector>::const_iterator data_set=cell_data.begin();
+ data_set!=cell_data.end(); ++data_set)
+ for (unsigned int i=0; i<data_set->names.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 <int dim>
+template <typename number>
+void DataOutStack<dim>::add_data_vector (const Vector<number> &vec,
+ const string &name)
+{
+ vector<string> names;
+ names.push_back (name);
+ add_data_vector (vec, names);
+};
+
+
+
+
+template <int dim>
+template <typename number>
+void DataOutStack<dim>::add_data_vector (const Vector<number> &vec,
+ const vector<string> &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<names.size(); ++i)
+ Assert (names[i].find_first_not_of("abcdefghijklmnopqrstuvwxyz"
+ "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+ "0123456789_<>()") == string::npos,
+ ExcInvalidCharacter (names[i]));
+
+ if (vec.size() == dof_handler->n_dofs())
+ {
+ vector<DataVector>::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<DataVector>::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 <int dim>
+void DataOutStack<dim>::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<dim>::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<dim> patch_points (q_trapez, n_subdivisions);
+ FEValues<dim> fe_patch_values (dof_handler->get_fe(),
+ patch_points,
+ update_default);
+ const unsigned int n_q_points = patch_points.n_quadrature_points;
+ vector<double> patch_values (n_q_points);
+ vector<Vector<double> > patch_values_system (n_q_points,
+ Vector<double>(n_components));
+
+ // add the required number of patches
+ DataOutBase::Patch<dim+1> 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<DataOutBase::Patch<dim+1> >::iterator patch = &patches[patches.size()-n_patches];
+ unsigned int cell_number = 0;
+ for (DoFHandler<dim>::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<dim+1>(cell->vertex(0)(0),
+ parameter-parameter_step);
+ patch->vertices[3] = Point<dim+1>(cell->vertex(0)(0),
+ parameter);
+ patch->vertices[1] = Point<dim+1>(cell->vertex(1)(0),
+ parameter-parameter_step);
+ patch->vertices[2] = Point<dim+1>(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; dataset<dof_data.size(); ++dataset)
+ {
+ if (n_components == 1)
+ {
+ fe_patch_values.get_function_values (dof_data[dataset].data,
+ patch_values);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ patch->data(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; component<n_components; ++component)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ patch->data(dataset*n_components+component,
+ q*(n_subdivisions+1)+i)
+ = patch_values_system[q](component);
+ };
+ };
+
+ // then do the cell data
+ for (unsigned int dataset=0; dataset<cell_data.size(); ++dataset)
+ {
+ const double value = cell_data[dataset].data(cell_number);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ patch->data(dataset+dof_data.size()*n_components,
+ q*(n_subdivisions+1)+i) = value;
+ };
+ };
+ };
+};
+
+
+
+template <int dim>
+void DataOutStack<dim>::finish_parameter_value ()
+{
+ // release lock on dof handler
+ dof_handler = 0;
+ for (vector<DataVector>::iterator i=dof_data.begin(); i!=dof_data.end(); ++i)
+ i->data.reinit (0);
+
+ for (vector<DataVector>::iterator i=cell_data.begin(); i!=cell_data.end(); ++i)
+ i->data.reinit (0);
+};
+
+
+
+template <int dim>
+const vector<DataOutBase::Patch<dim+1> > &
+DataOutStack<dim>::get_patches () const
+{
+ return patches;
+};
+
+
+
+template <int dim>
+vector<string> DataOutStack<dim>::get_dataset_names () const
+{
+ vector<string> names;
+ for (vector<DataVector>::const_iterator dataset=dof_data.begin();
+ dataset!=dof_data.end(); ++dataset)
+ names.insert (names.end(), dataset->names.begin(), dataset->names.end());
+ for (vector<DataVector>::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<deal_II_dimension>;
+template void DataOutStack<deal_II_dimension>::add_data_vector (const Vector<double> &vec,
+ const string &name);
+template void DataOutStack<deal_II_dimension>::add_data_vector (const Vector<float> &vec,
+ const string &name);
+