#include <base/data_out_base.h>
#include <dofs/dof_handler.h>
#include <fe/mapping.h>
+#include <numerics/data_postprocessor.h>
#include <boost/shared_ptr.hpp>
const std::string &name,
const DataVectorType type = type_automatic);
+ /**
+ * This function is an alternative to the
+ * above ones, allowing the output of
+ * derived quantities instead of the given
+ * data. This converison has to be done in
+ * a class derived from DataPostprocessor.
+ *
+ * The names for these derived quantities
+ * are implicitly given via the @p
+ * data_postprocessor argument. As only
+ * data of type @p type_dof_data can be
+ * transformed, this type is also known
+ * implicitly and does not have to be
+ * given.
+ *
+ * The actual type for the template
+ * argument may be any vector type from
+ * which FEValues can extract values
+ * on a cell using the
+ * FEValuesBase::get_function_values() function.
+ */
+ template <class VECTOR>
+ void add_data_vector (const VECTOR &data,
+ const DataPostprocessor<DH::dimension> &data_postprocessor);
+
/**
* Release the pointers to the
* data vectors. This allows
/**
* Exception
*/
+ DeclException0 (ExcDataPostprocessingIsNotPossibleForCellData);
+ /**
+ * Exception
+ */
DeclException3 (ExcInvalidVectorSize,
int, int, int,
<< "The vector has size " << arg1
* for the individual components of
* the vector.
*/
- DataEntryBase (const std::vector<std::string> &names);
+ DataEntryBase (const std::vector<std::string> &names,
+ const DataPostprocessor<DH::dimension> *data_postprocessor=0);
/**
* Destructor made virtual.
get_function_values (const FEValuesBase<DH::dimension> &fe_patch_values,
std::vector<Vector<double> > &patch_values_system) const = 0;
+ /**
+ * Given a FEValuesBase object,
+ * extract the gradients on the present
+ * cell from the vector we actually
+ * store.
+ */
+ virtual
+ void
+ get_function_gradients (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<Tensor<1,DH::dimension> > &patch_gradients) const = 0;
+
+ /**
+ * Given a FEValuesBase object,
+ * extract the gradients on the present
+ * cell from the vector we actually
+ * store. This function does the same
+ * as the one above but for
+ * vector-valued finite elements.
+ */
+ virtual
+ void
+ get_function_gradients (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<std::vector<Tensor<1,DH::dimension> > > &patch_gradients_system) const = 0;
+
+ /**
+ * Given a FEValuesBase object, extract
+ * the second derivatives on the
+ * present cell from the vector we
+ * actually store.
+ */
+ virtual
+ void
+ get_function_second_derivatives (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<Tensor<2,DH::dimension> > &patch_second_derivatives) const = 0;
+
+ /**
+ * Given a FEValuesBase object, extract
+ * the second derivatives on the
+ * present cell from the vector we
+ * actually store. This function does
+ * the same as the one above but for
+ * vector-valued finite elements.
+ */
+ virtual
+ void
+ get_function_second_derivatives (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<std::vector< Tensor<2,DH::dimension> > > &patch_second_derivatives_system) const = 0;
+
/**
* Clear all references to the
* vectors.
* Names of the components of this
* data vector.
*/
- std::vector<std::string> names;
+ const std::vector<std::string> names;
+
+ /**
+ * Pointer to a DataPostprocessing
+ * object which shall be applied to
+ * this data vector.
+ */
+ SmartPointer<const DataPostprocessor<DH::dimension> > postprocessor;
+
+ /**
+ * Number of output variables this
+ * dataset provides (either number of
+ * components in vector valued function
+ * / data vector or number of computed
+ * quantities, if DataPostprocessor is
+ * applied). This variable is
+ * determined via and thus equivalent
+ * to <tt>names.size()</tt>.
+ */
+ unsigned int n_output_variables;
};
* vector and the list of names for
* the individual components.
*/
- DataEntry (const VectorType *data,
- const std::vector<std::string> &names);
+ DataEntry (const VectorType *data,
+ const std::vector<std::string> &names,
+ const DataPostprocessor<DH::dimension> *data_postprocessor);
/**
* Assuming that the stored vector is
get_function_values (const FEValuesBase<DH::dimension> &fe_patch_values,
std::vector<Vector<double> > &patch_values_system) const;
+ /**
+ * Given a FEValuesBase object,
+ * extract the gradients on the present
+ * cell from the vector we actually
+ * store.
+ */
+ virtual
+ void
+ get_function_gradients (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<Tensor<1,DH::dimension> > &patch_gradients) const;
+
+ /**
+ * Given a FEValuesBase object,
+ * extract the gradients on the present
+ * cell from the vector we actually
+ * store. This function does the same
+ * as the one above but for
+ * vector-valued finite elements.
+ */
+ virtual
+ void
+ get_function_gradients (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<std::vector<Tensor<1,DH::dimension> > > &patch_gradients_system) const;
+
+ /**
+ * Given a FEValuesBase object, extract
+ * the second derivatives on the
+ * present cell from the vector we
+ * actually store.
+ */
+ virtual
+ void
+ get_function_second_derivatives (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<Tensor<2,DH::dimension> > &patch_second_derivatives) const;
+
+ /**
+ * Given a FEValuesBase object, extract
+ * the second derivatives on the
+ * present cell from the vector we
+ * actually store. This function does
+ * the same as the one above but for
+ * vector-valued finite elements.
+ */
+ virtual
+ void
+ get_function_second_derivatives (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<std::vector< Tensor<2,DH::dimension> > > &patch_second_derivatives_system) const;
+
/**
* Clear all references to the
* vectors.
*/
const VectorType *vector;
};
-
+
+ /**
+ * Common form of the three disctinct public
+ * functions with the same name, but
+ * different parameter list.
+ */
+ template <class VECTOR>
+ void add_data_vector (const VECTOR &data,
+ const std::vector<std::string> &names,
+ const DataVectorType type,
+ const DataPostprocessor<DH::dimension> *data_postprocessor);
+
+
/**
* Abbreviate the somewhat lengthy
SmartPointer<const Mapping<DH::dimension> > mapping;
std::vector<double> patch_values;
std::vector<Vector<double> > patch_values_system;
+ std::vector<Tensor<1,DH::dimension> > patch_gradients;
+ std::vector<std::vector<Tensor<1,DH::dimension> > > patch_gradients_system;
+ std::vector<Tensor<2,DH::dimension> > patch_second_derivatives;
+ std::vector<std::vector<Tensor<2,DH::dimension> > > patch_second_derivatives_system;
+ std::vector<Point<dim> > dummy_normals;
+ std::vector<std::vector<Vector<double> > > postprocessed_values;
std::vector<std::vector<unsigned int> > *cell_to_patch_index_map;
};
unsigned int n_components;
unsigned int n_datasets;
unsigned int n_subdivisions;
- std::vector<double> patch_values;
- std::vector<Vector<double> > patch_values_system;
+ std::vector<double> patch_values;
+ std::vector<Vector<double> > patch_values_system;
+ std::vector<Tensor<1,dim> > patch_gradients;
+ std::vector<std::vector<Tensor<1,dim> > > patch_gradients_system;
+ std::vector<Tensor<2,dim> > patch_second_derivatives;
+ std::vector<std::vector<Tensor<2,dim> > > patch_second_derivatives_system;
+ std::vector<std::vector<Vector<double> > > postprocessed_values;
+ std::vector<Point<dim> > patch_normals;
Data ()
{}
};
unsigned int n_datasets;
unsigned int n_patches_per_circle;
unsigned int n_subdivisions;
- std::vector<double> patch_values;
- std::vector<Vector<double> > patch_values_system;
+ std::vector<double> patch_values;
+ std::vector<Vector<double> > patch_values_system;
+ std::vector<Tensor<1,dim> > patch_gradients;
+ std::vector<std::vector<Tensor<1,dim> > > patch_gradients_system;
+ std::vector<Tensor<2,dim> > patch_second_derivatives;
+ std::vector<std::vector<Tensor<2,dim> > > patch_second_derivatives_system;
+ std::vector<std::vector<Vector<double> > > postprocessed_values;
+ std::vector<Point<dim> > dummy_normals;
Data ()
{}
};
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2007 by 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 __deal2__data_postprocessor_h
+#define __deal2__data_postprocessor_h
+
+
+
+#include <base/subscriptor.h>
+#include <base/tensor.h>
+#include <lac/vector.h>
+#include <fe/fe_update_flags.h>
+
+#include <vector>
+#include <string>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/*
+ * For the (graphical) output of a FE solution one frequently wants to include
+ * derived quantities, which are calculated from the values of the solution and
+ * possibly the first and second derivates of the solution. This class offers
+ * the interface to perform such a postprocessing. Given the values and
+ * derivatives of provided dataon given points of a cell, new quantities can be
+ * calculated.
+ *
+ * A data vector and an object of a derived class can be given to the
+ * <tt>DataOut::add_data_vector</tt> function, which will write the derived
+ * quantities instead of the provided data to the output file. Note, that the
+ * DataPostprocessor has to live until @p build_patches has been
+ * called. DataOutFaces and DataOutRotation can be used as well.
+ *
+ * In order not to perform needless calculations, DataPostprocessor has to
+ * provide the information, which input data is needed for the calculation of
+ * the derived quantities, i.e. whether it needs the values, the first
+ * derivative and/or the second derivative of the provided
+ * data. DataPostprocessor objects which are used in combination with a
+ * DataOutFaces object can also ask for the normal vectors at each point. The
+ * information, which data is needed has to be provided via the UpdateFlags
+ * returned by the virtual function @p get_needed_update_flags. It is your
+ * responsibility to use only those values which were updated in the calculation
+ * of derived quantities. The DataOut object will provide references to the
+ * requested data in the call to @p compute_derived_quantities.
+ *
+ * Furthermore, derived classes have to implement the @p get_names and @p
+ * n_output_variables functions, where the number of output variables returned
+ * by the latter function has to match the size of the vector returned by the
+ * former. Furthermore, this number has to match the number of computed
+ * quantities, of course.
+ *
+ * @ingroup output
+ * @author Tobias Leicht, 2007
+ */
+template <int dim>
+class DataPostprocessor: public Subscriptor
+{
+ public:
+ /**
+ * Virtual desctructor for safety. Does not
+ * do anything so far.
+ */
+ virtual ~DataPostprocessor ();
+
+ /**
+ * This is the main function which actually
+ * performs the postprocessing. The first
+ * argument is a reference to the
+ * postprocessed data which has correct
+ * size already and must be filled by this
+ * function. @p uh is a reference to a
+ * vector of data values at all points, @p
+ * duh the same for gradients, @p dduh for
+ * second derivatives and @p normals is a
+ * reference to the normal vectors. Note,
+ * that the last four references will only
+ * contain valid data, if the respective
+ * flags are returned by @p
+ * get_needed_update_flags, otherwise those
+ * vectors will be in an unspecified
+ * state. @p normals will always be an
+ * empty vector when working on cells, not
+ * on faces.
+ *
+ * This function is called, when the
+ * original data vector is scalar valued.
+ */
+ virtual void compute_derived_quantities (std::vector<Vector<double> > &computed_quantities,
+ const std::vector<double> &uh,
+ const std::vector<Tensor<1,dim> > &duh,
+ const std::vector<Tensor<2,dim> > &dduh,
+ const std::vector<Point<dim> > &normals) const;
+ /**
+ * Same as above function, but this one is
+ * called in case of vector-valued original
+ * data.
+ */
+ virtual void compute_derived_quantities (std::vector<Vector<double> > &computed_quantities,
+ const std::vector<Vector<double> > &uh,
+ const std::vector<std::vector<Tensor<1,dim> > > &duh,
+ const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+ const std::vector<Point<dim> > &normals) const;
+
+ /**
+ * Return the vector of strings describing
+ * the names of the computed quantities.
+ */
+ virtual std::vector<std::string> get_names () const=0;
+
+ /**
+ * Return, which data has to be provided to
+ * compute the derived quantities. This has
+ * to be a combination of @p update_values,
+ * @p update_gradients and @p
+ * update_second_derivatives. If the
+ * DataPostprocessor is to be used in
+ * combination with DataOutFaces, you may
+ * also ask for a update of normals via the
+ * @p update_normal_vectors flag.
+ */
+ virtual UpdateFlags get_needed_update_flags () const=0;
+
+ /**
+ * Number of postprocessed variables. Has
+ * to match the number of entries filled by
+ * @p compute_derived_quantities as well as
+ * the size of the vector of names returned
+ * by @p get_names.
+ */
+ virtual unsigned int n_output_variables() const=0;
+
+};
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
template <class DH, int patch_dim, int patch_space_dim>
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
-DataEntryBase::DataEntryBase (const std::vector<std::string> &names)
+DataEntryBase::DataEntryBase (const std::vector<std::string> &names_in,
+ const DataPostprocessor<DH::dimension> *data_postprocessor)
:
- names(names)
-{}
+ names(names_in),
+ postprocessor(data_postprocessor),
+ n_output_variables(names.size())
+{
+ Assert(!data_postprocessor ||
+ data_postprocessor->n_output_variables()==names.size(),
+ ExcDimensionMismatch(data_postprocessor->n_output_variables(),names.size()));
+}
template <typename VectorType>
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
DataEntry<VectorType>::
-DataEntry (const VectorType *data,
- const std::vector<std::string> &names)
+DataEntry (const VectorType *data,
+ const std::vector<std::string> &names,
+ const DataPostprocessor<DH::dimension> *data_postprocessor)
:
- DataEntryBase (names),
+ DataEntryBase (names,data_postprocessor),
vector (data)
{}
+template <class DH, int patch_dim, int patch_space_dim>
+template <typename VectorType>
+void
+DataOut_DoFData<DH,patch_dim,patch_space_dim>::
+DataEntry<VectorType>::
+get_function_gradients (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<std::vector<Tensor<1,DH::dimension> > > &patch_gradients_system) const
+{
+ fe_patch_values.get_function_grads (*vector, patch_gradients_system);
+}
+
+
+
+template <class DH,
+ int patch_dim, int patch_space_dim>
+template <typename VectorType>
+void
+DataOut_DoFData<DH,patch_dim,patch_space_dim>::
+DataEntry<VectorType>::
+get_function_gradients (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<Tensor<1,DH::dimension> > &patch_gradients) const
+{
+ fe_patch_values.get_function_grads (*vector, patch_gradients);
+}
+
+
+
+template <class DH, int patch_dim, int patch_space_dim>
+template <typename VectorType>
+void
+DataOut_DoFData<DH,patch_dim,patch_space_dim>::
+DataEntry<VectorType>::
+get_function_second_derivatives (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<std::vector<Tensor<2,DH::dimension> > > &patch_second_derivatives_system) const
+{
+ fe_patch_values.get_function_2nd_derivatives (*vector, patch_second_derivatives_system);
+}
+
+
+
+template <class DH,
+ int patch_dim, int patch_space_dim>
+template <typename VectorType>
+void
+DataOut_DoFData<DH,patch_dim,patch_space_dim>::
+DataEntry<VectorType>::
+get_function_second_derivatives (const FEValuesBase<DH::dimension> &fe_patch_values,
+ std::vector<Tensor<2,DH::dimension> > &patch_second_derivatives) const
+{
+ fe_patch_values.get_function_2nd_derivatives (*vector, patch_second_derivatives);
+}
+
+
+
template <class DH,
int patch_dim, int patch_space_dim>
template <typename VectorType>
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
-add_data_vector (const VECTOR &vec,
- const std::vector<std::string> &names,
- const DataVectorType type)
+add_data_vector (const VECTOR &vec,
+ const std::vector<std::string> &names,
+ const DataVectorType type,
+ const DataPostprocessor<DH::dimension> *data_postprocessor)
{
Assert (dofs != 0, ExcNoDoFHandlerSelected ());
ExcInvalidVectorSize (vec.size(),
dofs->n_dofs(),
dofs->get_tria().n_active_cells()));
- Assert (names.size() == dofs->get_fe().n_components(),
+ Assert (data_postprocessor || (names.size() == dofs->get_fe().n_components()),
ExcInvalidNumberOfNames (names.size(), dofs->get_fe().n_components()));
break;
// been handled above...
Assert (false, ExcInternalError());
};
-
+
+ // check names for invalid characters
for (unsigned int i=0; i<names.size(); ++i)
Assert (names[i].find_first_not_of("abcdefghijklmnopqrstuvwxyz"
"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
"0123456789_<>()")));
- DataEntryBase * new_entry = new DataEntry<VECTOR>(&vec, names);
+ DataEntryBase * new_entry = new DataEntry<VECTOR>(&vec, names, data_postprocessor);
if (actual_type == type_dof_data)
dof_data.push_back (boost::shared_ptr<DataEntryBase>(new_entry));
else
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
-add_data_vector (const VECTOR &vec,
- const std::string &name,
- const DataVectorType type)
+add_data_vector (const VECTOR &vec,
+ const std::string &name,
+ const DataVectorType type)
{
Assert (dofs != 0, ExcNoDoFHandlerSelected ());
const unsigned int n_components = dofs->get_fe().n_components ();
}
}
- add_data_vector (vec, names, type);
+ add_data_vector (vec, names, type, 0);
+}
+
+
+
+template <class DH,
+ int patch_dim, int patch_space_dim>
+template <class VECTOR>
+void
+DataOut_DoFData<DH,patch_dim,patch_space_dim>::
+add_data_vector (const VECTOR &vec,
+ const std::vector<std::string> &names,
+ const DataVectorType type)
+{
+ add_data_vector(vec,names,type,0);
+}
+
+
+
+template <class DH,
+ int patch_dim, int patch_space_dim>
+template <class VECTOR>
+void
+DataOut_DoFData<DH,patch_dim,patch_space_dim>::
+add_data_vector (const VECTOR &vec,
+ const DataPostprocessor<DH::dimension> &data_postprocessor)
+{
+ Assert (dofs != 0, ExcNoDoFHandlerSelected ());
+
+ std::vector<std::string> names=data_postprocessor.get_names();
+
+ add_data_vector (vec, names, type_dof_data, &data_postprocessor);
}
const hp::QCollection<DH::dimension> q_collection (patch_points);
const hp::FECollection<DH::dimension> fe_collection(this->dofs->get_fe());
const hp::MappingCollection<DH::dimension> mapping_collection(*(data.mapping));
+
+ UpdateFlags update_flags=update_values | update_q_points;
+ for (unsigned int i=0; i<this->dof_data.size(); ++i)
+ if (this->dof_data[i]->postprocessor)
+ update_flags |= this->dof_data[i]->postprocessor->get_needed_update_flags();
+ // perhaps update_normal_vectors is present,
+ // which would only be useful on faces, but
+ // we may not use it here.
+ Assert (!(update_flags & update_normal_vectors),
+ ExcMessage("The update of normal vectors may not be requested for evaluation of data on cells via DataPostprocessor."));
hp::FEValues<DH::dimension> x_fe_patch_values (mapping_collection,
fe_collection,
q_collection,
- update_values | update_q_points);
+ update_flags);
const unsigned int n_q_points = patch_points.n_quadrature_points;
patch->data(patch->data.size(0)-dim+i,q)=q_points[q][i];
}
+ // counter for data records
+ unsigned int offset=0;
+
// first fill dof_data
for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
{
- if (data.n_components == 1)
+ const DataPostprocessor<dim> *postprocessor=this->dof_data[dataset]->postprocessor;
+
+ if (postprocessor)
{
- this->dof_data[dataset]->get_function_values (fe_patch_values,
- data.patch_values);
-
+ // we have to postprocess the
+ // data, so determine, which
+ // fields have to be updated
+ const UpdateFlags update_flags=postprocessor->get_needed_update_flags();
+ if (data.n_components == 1)
+ {
+ // at each point there is
+ // only one component of
+ // value, gradient etc.
+ if (update_flags & update_values)
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values);
+ if (update_flags & update_gradients)
+ this->dof_data[dataset]->get_function_gradients (fe_patch_values,
+ data.patch_gradients);
+ if (update_flags & update_second_derivatives)
+ this->dof_data[dataset]->get_function_second_derivatives (fe_patch_values,
+ data.patch_second_derivatives);
+ postprocessor->compute_derived_quantities(data.postprocessed_values[dataset],
+ data.patch_values,
+ data.patch_gradients,
+ data.patch_second_derivatives,
+ data.dummy_normals);
+ }
+ else
+ {
+ // at each point there is
+ // a vector valued
+ // function and its
+ // derivative...
+ if (update_flags & update_values)
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values_system);
+ if (update_flags & update_gradients)
+ this->dof_data[dataset]->get_function_gradients (fe_patch_values,
+ data.patch_gradients_system);
+ if (update_flags & update_second_derivatives)
+ this->dof_data[dataset]->get_function_second_derivatives (fe_patch_values,
+ data.patch_second_derivatives_system);
+ postprocessor->compute_derived_quantities(data.postprocessed_values[dataset],
+ data.patch_values_system,
+ data.patch_gradients_system,
+ data.patch_second_derivatives_system,
+ data.dummy_normals);
+ }
+
for (unsigned int q=0; q<n_q_points; ++q)
- patch->data(dataset,q) = data.patch_values[q];
+ for (unsigned int component=0; component<this->dof_data[dataset]->n_output_variables;++component)
+ patch->data(offset+component,q)= data.postprocessed_values[dataset][q](component);
}
else
- // system of components
- {
- this->dof_data[dataset]->get_function_values (fe_patch_values,
- data.patch_values_system);
-
- for (unsigned int component=0; component<data.n_components;
- ++component)
- for (unsigned int q=0; q<n_q_points; ++q)
- patch->data(dataset*data.n_components+component,q) =
- data.patch_values_system[q](component);
- }
+ // now we use the given data
+ // vector without
+ // modifications. again, we
+ // treat single component
+ // functions separately for
+ // efficiency reasons.
+ if (data.n_components == 1)
+ {
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ patch->data(offset,q) = data.patch_values[q];
+ }
+ else
+ {
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values_system);
+ for (unsigned int component=0; component<data.n_components;
+ ++component)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ patch->data(offset+component,q) =
+ data.patch_values_system[q](component);
+ }
+ // increment the counter for the
+ // actual data record
+ offset+=this->dof_data[dataset]->n_output_variables;
}
// then do the cell data. only
const double value
= this->cell_data[dataset]->get_cell_data_value (cell_index);
for (unsigned int q=0; q<n_q_points; ++q)
- patch->data(dataset+this->dof_data.size()*data.n_components,q) =
+ patch->data(offset+dataset,q) =
value;
}
}
const unsigned int n_q_points = patch_points.n_quadrature_points;
const unsigned int n_components = this->dofs->get_fe().n_components();
- const unsigned int n_datasets = this->dof_data.size() * n_components +
- this->cell_data.size();
+ unsigned int n_datasets=this->cell_data.size();
+ for (unsigned int i=0; i<this->dof_data.size(); ++i)
+ n_datasets+= this->dof_data[i]->n_output_variables;
// clear the patches array
if (true)
thread_data[i].n_subdivisions = n_subdivisions;
thread_data[i].patch_values.resize (n_q_points);
thread_data[i].patch_values_system.resize (n_q_points);
+ thread_data[i].patch_gradients.resize (n_q_points);
+ thread_data[i].patch_gradients_system.resize (n_q_points);
+ thread_data[i].patch_second_derivatives.resize (n_q_points);
+ thread_data[i].patch_second_derivatives_system.resize (n_q_points);
+ thread_data[i].dummy_normals.clear();
+ thread_data[i].postprocessed_values.resize(this->dof_data.size());
thread_data[i].mapping = &mapping;
thread_data[i].cell_to_patch_index_map = &cell_to_patch_index_map;
for (unsigned int k=0; k<n_q_points; ++k)
- thread_data[i].patch_values_system[k].reinit(n_components);
+ {
+ thread_data[i].patch_values_system[k].reinit(n_components);
+ thread_data[i].patch_gradients_system[k].resize(n_components);
+ thread_data[i].patch_second_derivatives_system[k].resize(n_components);
+ }
+
+ for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
+ if (this->dof_data[dataset]->postprocessor)
+ thread_data[i].postprocessed_values[dataset].resize(n_q_points,Vector<double>(this->dof_data[dataset]->n_output_variables));
}
Threads::ThreadGroup<> threads;
DataOut_DoFData<DH,D2,D3>:: \
add_data_vector<V > (const V &, \
const std::vector<std::string> &, \
- const DataVectorType)
+ const DataVectorType); \
+\
+template void \
+DataOut_DoFData<DH,D2,D3>:: \
+add_data_vector<V > (const V &, \
+ const DataPostprocessor<DH::dimension> &)
#ifndef DEAL_II_USE_PETSC
# define INSTANTIATE_VECTORS(DH,D2,D3) \
const hp::QCollection<DH::dimension-1> q_collection (patch_points);
const hp::FECollection<DH::dimension> fe_collection(this->dofs->get_fe());
+ UpdateFlags update_flags=update_values;
+ for (unsigned int i=0; i<this->dof_data.size(); ++i)
+ if (this->dof_data[i]->postprocessor)
+ update_flags |= this->dof_data[i]->postprocessor->get_needed_update_flags();
+
hp::FEFaceValues<DH::dimension> x_fe_patch_values (fe_collection, q_collection,
- update_values);
+ update_flags);
const unsigned int n_q_points = patch_points.n_quadrature_points;
const FEFaceValues<DH::dimension> &fe_patch_values
= x_fe_patch_values.get_present_fe_values ();
+ // counter for data records
+ unsigned int offset=0;
+
// first fill dof_data
for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
{
- if (data.n_components == 1)
+ const DataPostprocessor<dim> *postprocessor=this->dof_data[dataset]->postprocessor;
+ if (postprocessor)
{
- this->dof_data[dataset]->get_function_values (fe_patch_values,
- data.patch_values);
-
+ // we have to postprocess the
+ // data, so determine, which
+ // fields have to be updated
+ const UpdateFlags update_flags=postprocessor->get_needed_update_flags();
+
+ // get normals, if
+ // needed. this is a
+ // geometrical information
+ // and thus does not depend
+ // on the number of
+ // components of the data
+ // vector
+ if (update_flags & update_normal_vectors)
+ data.patch_normals=fe_patch_values.get_normal_vectors();
+
+ if (data.n_components == 1)
+ {
+ // at each point there is
+ // only one component of
+ // value, gradient etc.
+ if (update_flags & update_values)
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values);
+ if (update_flags & update_gradients)
+ this->dof_data[dataset]->get_function_gradients (fe_patch_values,
+ data.patch_gradients);
+ if (update_flags & update_second_derivatives)
+ this->dof_data[dataset]->get_function_second_derivatives (fe_patch_values,
+ data.patch_second_derivatives);
+ postprocessor->compute_derived_quantities(data.postprocessed_values[dataset],
+ data.patch_values,
+ data.patch_gradients,
+ data.patch_second_derivatives,
+ data.patch_normals);
+ }
+ else
+ {
+ // at each point there is
+ // a vector valued
+ // function and its
+ // derivative...
+ if (update_flags & update_values)
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values_system);
+ if (update_flags & update_gradients)
+ this->dof_data[dataset]->get_function_gradients (fe_patch_values,
+ data.patch_gradients_system);
+ if (update_flags & update_second_derivatives)
+ this->dof_data[dataset]->get_function_second_derivatives (fe_patch_values,
+ data.patch_second_derivatives_system);
+ postprocessor->compute_derived_quantities(data.postprocessed_values[dataset],
+ data.patch_values_system,
+ data.patch_gradients_system,
+ data.patch_second_derivatives_system,
+ data.patch_normals);
+ }
+
for (unsigned int q=0; q<n_q_points; ++q)
- patch->data(dataset,q) = data.patch_values[q];
+ for (unsigned int component=0; component<this->dof_data[dataset]->n_output_variables;++component)
+ patch->data(offset+component,q)= data.postprocessed_values[dataset][q](component);
}
else
- // system of components
- {
- this->dof_data[dataset]->get_function_values (fe_patch_values,
- data.patch_values_system);
-
- for (unsigned int component=0; component<data.n_components;
- ++component)
- for (unsigned int q=0; q<n_q_points; ++q)
- patch->data(dataset*data.n_components+component,q) =
- data.patch_values_system[q](component);
- };
- };
+ // now we use the given data
+ // vector without
+ // modifications. again, we
+ // treat single component
+ // functions separately for
+ // efficiency reasons.
+ if (data.n_components == 1)
+ {
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ patch->data(offset,q) = data.patch_values[q];
+ }
+ else
+ {
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values_system);
+ for (unsigned int component=0; component<data.n_components;
+ ++component)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ patch->data(offset+component,q) =
+ data.patch_values_system[q](component);
+ }
+ // increment the counter for the
+ // actual data record
+ offset+=this->dof_data[dataset]->n_output_variables;
+ }
// then do the cell data
for (unsigned int dataset=0; dataset<this->cell_data.size(); ++dataset)
const double value
= this->cell_data[dataset]->get_cell_data_value (cell_number);
for (unsigned int q=0; q<n_q_points; ++q)
- patch->data(dataset+this->dof_data.size()*data.n_components,q) =
+ patch->data(dataset+offset,q) =
value;
- };
- };
+ }
+ }
// next cell (patch) in this thread
for (unsigned int i=0;
(i<data.n_threads)&&(face.first != this->dofs->end()); ++i)
++patch;
face=next_face(face);
}
- };
+ }
}
const unsigned int n_q_points = patch_points.n_quadrature_points;
const unsigned int n_components = this->dofs->get_fe().n_components();
- const unsigned int n_datasets = this->dof_data.size() * n_components +
- this->cell_data.size();
+ unsigned int n_datasets = this->cell_data.size();
+ for (unsigned int i=0; i<this->dof_data.size(); ++i)
+ n_datasets+= this->dof_data[i]->n_output_variables;
// clear the patches array
if (true)
thread_data[i].n_subdivisions = n_subdivisions;
thread_data[i].patch_values.resize (n_q_points);
thread_data[i].patch_values_system.resize (n_q_points);
+ thread_data[i].patch_gradients.resize (n_q_points);
+ thread_data[i].patch_gradients_system.resize (n_q_points);
+ thread_data[i].patch_second_derivatives.resize (n_q_points);
+ thread_data[i].patch_second_derivatives_system.resize (n_q_points);
+ thread_data[i].patch_normals.resize (n_q_points);
+ thread_data[i].postprocessed_values.resize (this->dof_data.size());
for (unsigned int k=0; k<n_q_points; ++k)
- thread_data[i].patch_values_system[k].reinit(n_components);
+ {
+ thread_data[i].patch_values_system[k].reinit(n_components);
+ thread_data[i].patch_gradients_system[k].resize(n_components);
+ thread_data[i].patch_second_derivatives_system[k].resize(n_components);
+ }
+
+ for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
+ if (this->dof_data[dataset]->postprocessor)
+ thread_data[i].postprocessed_values[dataset].resize(n_q_points,Vector<double>(this->dof_data[dataset]->n_output_variables));
}
// create the patches with default
const hp::QCollection<DH::dimension> q_collection (patch_points);
const hp::FECollection<DH::dimension> fe_collection(this->dofs->get_fe());
+ UpdateFlags update_flags=update_values | update_q_points;
+ for (unsigned int i=0; i<this->dof_data.size(); ++i)
+ if (this->dof_data[i]->postprocessor)
+ update_flags |= this->dof_data[i]->postprocessor->get_needed_update_flags();
+ // perhaps update_normal_vectors is present,
+ // which would only be useful on faces, but
+ // we may not use it here.
+ Assert (!(update_flags & update_normal_vectors),
+ ExcMessage("The update of normal vectors may not be requested for evaluation of data on cells via DataPostprocessor."));
+
hp::FEValues<DH::dimension> x_fe_patch_values (fe_collection, q_collection,
- update_values);
+ update_flags);
const unsigned int n_patches_per_circle = data.n_patches_per_circle;
Assert (false, ExcNotImplemented());
};
+ unsigned int offset=0;
// then fill in data
if (data.n_datasets > 0)
// first fill dof_data
for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
{
- if (data.n_components == 1)
+ const DataPostprocessor<dim> *postprocessor=this->dof_data[dataset]->postprocessor;
+ if (postprocessor)
{
- this->dof_data[dataset]->get_function_values (fe_patch_values,
- data.patch_values);
-
- switch (DH::dimension)
+ // we have to postprocess the
+ // data, so determine, which
+ // fields have to be updated
+ const UpdateFlags update_flags=postprocessor->get_needed_update_flags();
+
+ if (data.n_components == 1)
{
- case 1:
- for (unsigned int x=0; x<n_points; ++x)
- for (unsigned int y=0; y<n_points; ++y)
- patch->data(dataset,
- x*n_points + y)
- = data.patch_values[x];
- break;
-
- case 2:
- for (unsigned int x=0; x<n_points; ++x)
- for (unsigned int y=0; y<n_points; ++y)
- for (unsigned int z=0; z<n_points; ++z)
- patch->data(dataset,
- x*n_points*n_points +
- y +
- z*n_points)
- = data.patch_values[x*n_points+z];
- break;
-
- default:
- Assert (false, ExcNotImplemented());
- };
+ // at each point there is
+ // only one component of
+ // value, gradient etc.
+ if (update_flags & update_values)
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values);
+ if (update_flags & update_gradients)
+ this->dof_data[dataset]->get_function_gradients (fe_patch_values,
+ data.patch_gradients);
+ if (update_flags & update_second_derivatives)
+ this->dof_data[dataset]->get_function_second_derivatives (fe_patch_values,
+ data.patch_second_derivatives);
+ postprocessor->compute_derived_quantities(data.postprocessed_values[dataset],
+ data.patch_values,
+ data.patch_gradients,
+ data.patch_second_derivatives,
+ data.dummy_normals);
+ }
+ else
+ {
+ // at each point there is
+ // a vector valued
+ // function and its
+ // derivative...
+ if (update_flags & update_values)
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values_system);
+ if (update_flags & update_gradients)
+ this->dof_data[dataset]->get_function_gradients (fe_patch_values,
+ data.patch_gradients_system);
+ if (update_flags & update_second_derivatives)
+ this->dof_data[dataset]->get_function_second_derivatives (fe_patch_values,
+ data.patch_second_derivatives_system);
+ postprocessor->compute_derived_quantities(data.postprocessed_values[dataset],
+ data.patch_values_system,
+ data.patch_gradients_system,
+ data.patch_second_derivatives_system,
+ data.dummy_normals);
+ }
+
+ for (unsigned int component=0;
+ component<this->dof_data[dataset]->n_output_variables;
+ ++component)
+ {
+ switch (DH::dimension)
+ {
+ case 1:
+ for (unsigned int x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ patch->data(offset+component,
+ x*n_points + y)
+ = data.postprocessed_values[dataset][x](component);
+ break;
+
+ case 2:
+ for (unsigned int x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ for (unsigned int z=0; z<n_points; ++z)
+ patch->data(offset+component,
+ x*n_points*n_points +
+ y*n_points +
+ z)
+ = data.postprocessed_values[dataset][x*n_points+z](component);
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
}
else
- // system of components
- {
- this->dof_data[dataset]->get_function_values (fe_patch_values,
- data.patch_values_system);
-
- for (unsigned int component=0; component<data.n_components;
- ++component)
- {
- switch (DH::dimension)
- {
- case 1:
- for (unsigned int x=0; x<n_points; ++x)
- for (unsigned int y=0; y<n_points; ++y)
- patch->data(dataset*data.n_components+component,
- x*n_points + y)
+ if (data.n_components == 1)
+ {
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values);
+
+ switch (DH::dimension)
+ {
+ case 1:
+ for (unsigned int x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ patch->data(offset,
+ x*n_points + y)
+ = data.patch_values[x];
+ break;
+
+ case 2:
+ for (unsigned int x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ for (unsigned int z=0; z<n_points; ++z)
+ patch->data(offset,
+ x*n_points*n_points +
+ y +
+ z*n_points)
+ = data.patch_values[x*n_points+z];
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ else
+ // system of components
+ {
+ this->dof_data[dataset]->get_function_values (fe_patch_values,
+ data.patch_values_system);
+
+ for (unsigned int component=0; component<data.n_components;
+ ++component)
+ {
+ switch (DH::dimension)
+ {
+ case 1:
+ for (unsigned int x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ patch->data(offset+component,
+ x*n_points + y)
= data.patch_values_system[x](component);
- break;
-
- case 2:
- for (unsigned int x=0; x<n_points; ++x)
- for (unsigned int y=0; y<n_points; ++y)
- for (unsigned int z=0; z<n_points; ++z)
- patch->data(dataset*data.n_components+component,
- x*n_points*n_points +
- y*n_points +
- z)
- = data.patch_values_system[x*n_points+z](component);
- break;
-
- default:
- Assert (false, ExcNotImplemented());
- };
- };
- };
- };
+ break;
+
+ case 2:
+ for (unsigned int x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ for (unsigned int z=0; z<n_points; ++z)
+ patch->data(offset+component,
+ x*n_points*n_points +
+ y*n_points +
+ z)
+ = data.patch_values_system[x*n_points+z](component);
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ }
+ offset+=this->dof_data[dataset]->n_output_variables;
+ }
// then do the cell data
for (unsigned int dataset=0; dataset<this->cell_data.size(); ++dataset)
case 1:
for (unsigned int x=0; x<n_points; ++x)
for (unsigned int y=0; y<n_points; ++y)
- patch->data(dataset+this->dof_data.size()*data.n_components,
+ patch->data(dataset+offset,
x*n_points +
y)
= value;
for (unsigned int x=0; x<n_points; ++x)
for (unsigned int y=0; y<n_points; ++y)
for (unsigned int z=0; z<n_points; ++z)
- patch->data(dataset+this->dof_data.size()*data.n_components,
+ patch->data(dataset+offset,
x*n_points*n_points +
y*n_points +
z)
default:
Assert (false, ExcNotImplemented());
- };
- };
- };
- };
+ }
+ }
+ }
+ }
// next cell (patch) in this
// thread. note that we have
{
++cell_number;
cell=next_cell(cell);
- };
- };
+ }
+ }
}
const unsigned int n_q_points = patch_points.n_quadrature_points;
const unsigned int n_components = this->dofs->get_fe().n_components();
- const unsigned int n_datasets = this->dof_data.size() * n_components +
- this->cell_data.size();
+ unsigned int n_datasets=this->cell_data.size();
+ for (unsigned int i=0; i<this->dof_data.size(); ++i)
+ n_datasets+= this->dof_data[i]->n_output_variables;
// clear the patches array
if (true)
thread_data[i].n_subdivisions = n_subdivisions;
thread_data[i].patch_values.resize (n_q_points);
thread_data[i].patch_values_system.resize (n_q_points);
+ thread_data[i].patch_gradients.resize (n_q_points);
+ thread_data[i].patch_gradients_system.resize (n_q_points);
+ thread_data[i].patch_second_derivatives.resize (n_q_points);
+ thread_data[i].patch_second_derivatives_system.resize (n_q_points);
+ thread_data[i].dummy_normals.clear();
+ thread_data[i].postprocessed_values.resize(this->dof_data.size());
for (unsigned int k=0; k<n_q_points; ++k)
- thread_data[i].patch_values_system[k].reinit(n_components);
+ {
+ thread_data[i].patch_values_system[k].reinit(n_components);
+ thread_data[i].patch_gradients_system[k].resize(n_components);
+ thread_data[i].patch_second_derivatives_system[k].resize(n_components);
+ }
+
+ for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
+ if (this->dof_data[dataset]->postprocessor)
+ thread_data[i].postprocessed_values[dataset].resize(n_q_points,Vector<double>(this->dof_data[dataset]->n_output_variables));
}
// create the patches with default
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2007 by 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 <numerics/data_postprocessor.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template <int dim>
+DataPostprocessor<dim>::~DataPostprocessor()
+{}
+
+
+
+template <int dim>
+void
+DataPostprocessor<dim>::compute_derived_quantities (std::vector<Vector<double> > &computed_quantities,
+ const std::vector<double> &/*uh*/,
+ const std::vector<Tensor<1,dim> > &/*duh*/,
+ const std::vector<Tensor<2,dim> > &/*dduh*/,
+ const std::vector<Point<dim> > &/*normals*/) const
+{
+ computed_quantities.clear();
+ AssertThrow(false,ExcPureFunctionCalled());
+}
+
+
+
+template <int dim>
+void
+DataPostprocessor<dim>::compute_derived_quantities (std::vector<Vector<double> > &computed_quantities,
+ const std::vector<Vector<double> > &/*uh*/,
+ const std::vector<std::vector<Tensor<1,dim> > > &/*duh*/,
+ const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+ const std::vector<Point<dim> > &/*normals*/) const
+{
+ computed_quantities.clear();
+ AssertThrow(false,ExcPureFunctionCalled());
+}
+
+// explicit instantiation
+template class DataPostprocessor<deal_II_dimension>;
+
+DEAL_II_NAMESPACE_CLOSE