From 6883c7b61f65208944cc813952027525fdbad7e8 Mon Sep 17 00:00:00 2001 From: leicht Date: Mon, 10 Sep 2007 11:22:30 +0000 Subject: [PATCH] Introduce new class DataPostprocessor, which takes data values and possibly gradients, second derivatives (and normals on faces) and calculates derived data from a given FE solution. Derived classes can be used in combination with DataOut, DataOutFaces and DataOutRotation. git-svn-id: https://svn.dealii.org/trunk@15186 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/data_out.h | 175 +++++++++++- .../deal.II/include/numerics/data_out_faces.h | 10 +- .../include/numerics/data_out_rotation.h | 10 +- .../include/numerics/data_postprocessor.h | 145 ++++++++++ deal.II/deal.II/source/numerics/data_out.cc | 265 +++++++++++++++--- .../deal.II/source/numerics/data_out_faces.cc | 144 ++++++++-- .../source/numerics/data_out_rotation.cc | 254 ++++++++++++----- .../source/numerics/data_postprocessor.cc | 53 ++++ 8 files changed, 911 insertions(+), 145 deletions(-) create mode 100644 deal.II/deal.II/include/numerics/data_postprocessor.h create mode 100644 deal.II/deal.II/source/numerics/data_postprocessor.cc diff --git a/deal.II/deal.II/include/numerics/data_out.h b/deal.II/deal.II/include/numerics/data_out.h index 502060100b..94cec56ec0 100644 --- a/deal.II/deal.II/include/numerics/data_out.h +++ b/deal.II/deal.II/include/numerics/data_out.h @@ -21,6 +21,7 @@ #include #include #include +#include #include @@ -343,6 +344,31 @@ class DataOut_DoFData : public DataOutInterface 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 + void add_data_vector (const VECTOR &data, + const DataPostprocessor &data_postprocessor); + /** * Release the pointers to the * data vectors. This allows @@ -460,6 +486,10 @@ class DataOut_DoFData : public DataOutInterface /** * Exception */ + DeclException0 (ExcDataPostprocessingIsNotPossibleForCellData); + /** + * Exception + */ DeclException3 (ExcInvalidVectorSize, int, int, int, << "The vector has size " << arg1 @@ -533,7 +563,8 @@ class DataOut_DoFData : public DataOutInterface * for the individual components of * the vector. */ - DataEntryBase (const std::vector &names); + DataEntryBase (const std::vector &names, + const DataPostprocessor *data_postprocessor=0); /** * Destructor made virtual. @@ -573,6 +604,54 @@ class DataOut_DoFData : public DataOutInterface get_function_values (const FEValuesBase &fe_patch_values, std::vector > &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 &fe_patch_values, + std::vector > &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 &fe_patch_values, + std::vector > > &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 &fe_patch_values, + std::vector > &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 &fe_patch_values, + std::vector > > &patch_second_derivatives_system) const = 0; + /** * Clear all references to the * vectors. @@ -590,7 +669,26 @@ class DataOut_DoFData : public DataOutInterface * Names of the components of this * data vector. */ - std::vector names; + const std::vector names; + + /** + * Pointer to a DataPostprocessing + * object which shall be applied to + * this data vector. + */ + SmartPointer > 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 names.size(). + */ + unsigned int n_output_variables; }; @@ -611,8 +709,9 @@ class DataOut_DoFData : public DataOutInterface * vector and the list of names for * the individual components. */ - DataEntry (const VectorType *data, - const std::vector &names); + DataEntry (const VectorType *data, + const std::vector &names, + const DataPostprocessor *data_postprocessor); /** * Assuming that the stored vector is @@ -647,6 +746,54 @@ class DataOut_DoFData : public DataOutInterface get_function_values (const FEValuesBase &fe_patch_values, std::vector > &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 &fe_patch_values, + std::vector > &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 &fe_patch_values, + std::vector > > &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 &fe_patch_values, + std::vector > &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 &fe_patch_values, + std::vector > > &patch_second_derivatives_system) const; + /** * Clear all references to the * vectors. @@ -670,7 +817,19 @@ class DataOut_DoFData : public DataOutInterface */ const VectorType *vector; }; - + + /** + * Common form of the three disctinct public + * functions with the same name, but + * different parameter list. + */ + template + void add_data_vector (const VECTOR &data, + const std::vector &names, + const DataVectorType type, + const DataPostprocessor *data_postprocessor); + + /** * Abbreviate the somewhat lengthy @@ -982,6 +1141,12 @@ class DataOut : public DataOut_DoFData SmartPointer > mapping; std::vector patch_values; std::vector > patch_values_system; + std::vector > patch_gradients; + std::vector > > patch_gradients_system; + std::vector > patch_second_derivatives; + std::vector > > patch_second_derivatives_system; + std::vector > dummy_normals; + std::vector > > postprocessed_values; std::vector > *cell_to_patch_index_map; }; diff --git a/deal.II/deal.II/include/numerics/data_out_faces.h b/deal.II/deal.II/include/numerics/data_out_faces.h index ca124d9d93..49b17a46d2 100644 --- a/deal.II/deal.II/include/numerics/data_out_faces.h +++ b/deal.II/deal.II/include/numerics/data_out_faces.h @@ -196,8 +196,14 @@ class DataOutFaces : public DataOut_DoFData patch_values; - std::vector > patch_values_system; + std::vector patch_values; + std::vector > patch_values_system; + std::vector > patch_gradients; + std::vector > > patch_gradients_system; + std::vector > patch_second_derivatives; + std::vector > > patch_second_derivatives_system; + std::vector > > postprocessed_values; + std::vector > patch_normals; Data () {} }; diff --git a/deal.II/deal.II/include/numerics/data_out_rotation.h b/deal.II/deal.II/include/numerics/data_out_rotation.h index b02ffdf0e0..c8733d1228 100644 --- a/deal.II/deal.II/include/numerics/data_out_rotation.h +++ b/deal.II/deal.II/include/numerics/data_out_rotation.h @@ -185,8 +185,14 @@ class DataOutRotation : public DataOut_DoFData unsigned int n_datasets; unsigned int n_patches_per_circle; unsigned int n_subdivisions; - std::vector patch_values; - std::vector > patch_values_system; + std::vector patch_values; + std::vector > patch_values_system; + std::vector > patch_gradients; + std::vector > > patch_gradients_system; + std::vector > patch_second_derivatives; + std::vector > > patch_second_derivatives_system; + std::vector > > postprocessed_values; + std::vector > dummy_normals; Data () {} }; diff --git a/deal.II/deal.II/include/numerics/data_postprocessor.h b/deal.II/deal.II/include/numerics/data_postprocessor.h new file mode 100644 index 0000000000..553fec02bb --- /dev/null +++ b/deal.II/deal.II/include/numerics/data_postprocessor.h @@ -0,0 +1,145 @@ +//--------------------------------------------------------------------------- +// $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 +#include +#include +#include + +#include +#include + +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 + * DataOut::add_data_vector 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 +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 > &computed_quantities, + const std::vector &uh, + const std::vector > &duh, + const std::vector > &dduh, + const std::vector > &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 > &computed_quantities, + const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals) const; + + /** + * Return the vector of strings describing + * the names of the computed quantities. + */ + virtual std::vector 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 diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index 56ab3c83a5..db12a05d03 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -37,10 +37,17 @@ DEAL_II_NAMESPACE_OPEN template DataOut_DoFData:: -DataEntryBase::DataEntryBase (const std::vector &names) +DataEntryBase::DataEntryBase (const std::vector &names_in, + const DataPostprocessor *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())); +} @@ -55,10 +62,11 @@ template template DataOut_DoFData:: DataEntry:: -DataEntry (const VectorType *data, - const std::vector &names) +DataEntry (const VectorType *data, + const std::vector &names, + const DataPostprocessor *data_postprocessor) : - DataEntryBase (names), + DataEntryBase (names,data_postprocessor), vector (data) {} @@ -102,6 +110,60 @@ get_function_values (const FEValuesBase &fe_patch_values, +template +template +void +DataOut_DoFData:: +DataEntry:: +get_function_gradients (const FEValuesBase &fe_patch_values, + std::vector > > &patch_gradients_system) const +{ + fe_patch_values.get_function_grads (*vector, patch_gradients_system); +} + + + +template +template +void +DataOut_DoFData:: +DataEntry:: +get_function_gradients (const FEValuesBase &fe_patch_values, + std::vector > &patch_gradients) const +{ + fe_patch_values.get_function_grads (*vector, patch_gradients); +} + + + +template +template +void +DataOut_DoFData:: +DataEntry:: +get_function_second_derivatives (const FEValuesBase &fe_patch_values, + std::vector > > &patch_second_derivatives_system) const +{ + fe_patch_values.get_function_2nd_derivatives (*vector, patch_second_derivatives_system); +} + + + +template +template +void +DataOut_DoFData:: +DataEntry:: +get_function_second_derivatives (const FEValuesBase &fe_patch_values, + std::vector > &patch_second_derivatives) const +{ + fe_patch_values.get_function_2nd_derivatives (*vector, patch_second_derivatives); +} + + + template template @@ -162,9 +224,10 @@ template void DataOut_DoFData:: -add_data_vector (const VECTOR &vec, - const std::vector &names, - const DataVectorType type) +add_data_vector (const VECTOR &vec, + const std::vector &names, + const DataVectorType type, + const DataPostprocessor *data_postprocessor) { Assert (dofs != 0, ExcNoDoFHandlerSelected ()); @@ -203,7 +266,7 @@ add_data_vector (const VECTOR &vec, 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; @@ -212,7 +275,8 @@ add_data_vector (const VECTOR &vec, // been handled above... Assert (false, ExcInternalError()); }; - + + // check names for invalid characters for (unsigned int i=0; i()"))); - DataEntryBase * new_entry = new DataEntry(&vec, names); + DataEntryBase * new_entry = new DataEntry(&vec, names, data_postprocessor); if (actual_type == type_dof_data) dof_data.push_back (boost::shared_ptr(new_entry)); else @@ -236,9 +300,9 @@ template void DataOut_DoFData:: -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 (); @@ -265,7 +329,38 @@ add_data_vector (const VECTOR &vec, } } - add_data_vector (vec, names, type); + add_data_vector (vec, names, type, 0); +} + + + +template +template +void +DataOut_DoFData:: +add_data_vector (const VECTOR &vec, + const std::vector &names, + const DataVectorType type) +{ + add_data_vector(vec,names,type,0); +} + + + +template +template +void +DataOut_DoFData:: +add_data_vector (const VECTOR &vec, + const DataPostprocessor &data_postprocessor) +{ + Assert (dofs != 0, ExcNoDoFHandlerSelected ()); + + std::vector names=data_postprocessor.get_names(); + + add_data_vector (vec, names, type_dof_data, &data_postprocessor); } @@ -401,11 +496,21 @@ void DataOut::build_some_patches (Data &data) const hp::QCollection q_collection (patch_points); const hp::FECollection fe_collection(this->dofs->get_fe()); const hp::MappingCollection mapping_collection(*(data.mapping)); + + UpdateFlags update_flags=update_values | update_q_points; + for (unsigned int i=0; idof_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 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; @@ -483,29 +588,93 @@ void DataOut::build_some_patches (Data &data) 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; datasetdof_data.size(); ++dataset) { - if (data.n_components == 1) + const DataPostprocessor *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; qdata(dataset,q) = data.patch_values[q]; + for (unsigned int component=0; componentdof_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; componentdata(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; qdata(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; componentdata(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 @@ -523,7 +692,7 @@ void DataOut::build_some_patches (Data &data) const double value = this->cell_data[dataset]->get_cell_data_value (cell_index); for (unsigned int q=0; qdata(dataset+this->dof_data.size()*data.n_components,q) = + patch->data(offset+dataset,q) = value; } } @@ -642,8 +811,9 @@ void DataOut::build_patches (const Mapping &mapping, 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; idof_data.size(); ++i) + n_datasets+= this->dof_data[i]->n_output_variables; // clear the patches array if (true) @@ -719,12 +889,26 @@ void DataOut::build_patches (const Mapping &mapping, 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; kdof_data.size(); ++dataset) + if (this->dof_data[dataset]->postprocessor) + thread_data[i].postprocessed_values[dataset].resize(n_q_points,Vector(this->dof_data[dataset]->n_output_variables)); } Threads::ThreadGroup<> threads; @@ -770,7 +954,12 @@ template void \ DataOut_DoFData:: \ add_data_vector (const V &, \ const std::vector &, \ - const DataVectorType) + const DataVectorType); \ +\ +template void \ +DataOut_DoFData:: \ +add_data_vector (const V &, \ + const DataPostprocessor &) #ifndef DEAL_II_USE_PETSC # define INSTANTIATE_VECTORS(DH,D2,D3) \ diff --git a/deal.II/deal.II/source/numerics/data_out_faces.cc b/deal.II/deal.II/source/numerics/data_out_faces.cc index bb0eed7fe1..20b3449f2b 100644 --- a/deal.II/deal.II/source/numerics/data_out_faces.cc +++ b/deal.II/deal.II/source/numerics/data_out_faces.cc @@ -58,8 +58,13 @@ void DataOutFaces::build_some_patches (Data &data) const hp::QCollection q_collection (patch_points); const hp::FECollection fe_collection(this->dofs->get_fe()); + UpdateFlags update_flags=update_values; + for (unsigned int i=0; idof_data.size(); ++i) + if (this->dof_data[i]->postprocessor) + update_flags |= this->dof_data[i]->postprocessor->get_needed_update_flags(); + hp::FEFaceValues x_fe_patch_values (fe_collection, q_collection, - update_values); + update_flags); const unsigned int n_q_points = patch_points.n_quadrature_points; @@ -88,30 +93,104 @@ void DataOutFaces::build_some_patches (Data &data) const FEFaceValues &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; datasetdof_data.size(); ++dataset) { - if (data.n_components == 1) + const DataPostprocessor *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; qdata(dataset,q) = data.patch_values[q]; + for (unsigned int component=0; componentdof_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; componentdata(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; qdata(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; componentdata(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; datasetcell_data.size(); ++dataset) @@ -134,10 +213,10 @@ void DataOutFaces::build_some_patches (Data &data) const double value = this->cell_data[dataset]->get_cell_data_value (cell_number); for (unsigned int q=0; qdata(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; (idofs->end()); ++i) @@ -145,7 +224,7 @@ void DataOutFaces::build_some_patches (Data &data) ++patch; face=next_face(face); } - }; + } } @@ -178,8 +257,9 @@ void DataOutFaces::build_patches (const unsigned int nnnn_subdivisions, 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; idof_data.size(); ++i) + n_datasets+= this->dof_data[i]->n_output_variables; // clear the patches array if (true) @@ -209,9 +289,23 @@ void DataOutFaces::build_patches (const unsigned int nnnn_subdivisions, 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; kdof_data.size(); ++dataset) + if (this->dof_data[dataset]->postprocessor) + thread_data[i].postprocessed_values[dataset].resize(n_q_points,Vector(this->dof_data[dataset]->n_output_variables)); } // create the patches with default diff --git a/deal.II/deal.II/source/numerics/data_out_rotation.cc b/deal.II/deal.II/source/numerics/data_out_rotation.cc index 479e8002e3..6b2059a0d5 100644 --- a/deal.II/deal.II/source/numerics/data_out_rotation.cc +++ b/deal.II/deal.II/source/numerics/data_out_rotation.cc @@ -64,8 +64,18 @@ void DataOutRotation::build_some_patches (Data &data) const hp::QCollection q_collection (patch_points); const hp::FECollection fe_collection(this->dofs->get_fe()); + UpdateFlags update_flags=update_values | update_q_points; + for (unsigned int i=0; idof_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 x_fe_patch_values (fe_collection, q_collection, - update_values); + update_flags); const unsigned int n_patches_per_circle = data.n_patches_per_circle; @@ -171,6 +181,7 @@ void DataOutRotation::build_some_patches (Data &data) Assert (false, ExcNotImplemented()); }; + unsigned int offset=0; // then fill in data if (data.n_datasets > 0) @@ -182,72 +193,154 @@ void DataOutRotation::build_some_patches (Data &data) // first fill dof_data for (unsigned int dataset=0; datasetdof_data.size(); ++dataset) { - if (data.n_components == 1) + const DataPostprocessor *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; xdata(dataset, - x*n_points + y) - = data.patch_values[x]; - break; - - case 2: - for (unsigned int x=0; xdata(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; + componentdof_data[dataset]->n_output_variables; + ++component) + { + switch (DH::dimension) + { + case 1: + for (unsigned int x=0; xdata(offset+component, + x*n_points + y) + = data.postprocessed_values[dataset][x](component); + break; + + case 2: + for (unsigned int x=0; xdata(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; componentdata(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; xdata(offset, + x*n_points + y) + = data.patch_values[x]; + break; + + case 2: + for (unsigned int x=0; xdata(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; componentdata(offset+component, + x*n_points + y) = data.patch_values_system[x](component); - break; - - case 2: - for (unsigned int x=0; xdata(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; xdata(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; datasetcell_data.size(); ++dataset) @@ -259,7 +352,7 @@ void DataOutRotation::build_some_patches (Data &data) case 1: for (unsigned int x=0; xdata(dataset+this->dof_data.size()*data.n_components, + patch->data(dataset+offset, x*n_points + y) = value; @@ -269,7 +362,7 @@ void DataOutRotation::build_some_patches (Data &data) for (unsigned int x=0; xdata(dataset+this->dof_data.size()*data.n_components, + patch->data(dataset+offset, x*n_points*n_points + y*n_points + z) @@ -278,10 +371,10 @@ void DataOutRotation::build_some_patches (Data &data) default: Assert (false, ExcNotImplemented()); - }; - }; - }; - }; + } + } + } + } // next cell (patch) in this // thread. note that we have @@ -302,8 +395,8 @@ void DataOutRotation::build_some_patches (Data &data) { ++cell_number; cell=next_cell(cell); - }; - }; + } + } } @@ -357,8 +450,9 @@ void DataOutRotation::build_patches ( 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; idof_data.size(); ++i) + n_datasets+= this->dof_data[i]->n_output_variables; // clear the patches array if (true) @@ -394,9 +488,23 @@ void DataOutRotation::build_patches ( 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; kdof_data.size(); ++dataset) + if (this->dof_data[dataset]->postprocessor) + thread_data[i].postprocessed_values[dataset].resize(n_q_points,Vector(this->dof_data[dataset]->n_output_variables)); } // create the patches with default diff --git a/deal.II/deal.II/source/numerics/data_postprocessor.cc b/deal.II/deal.II/source/numerics/data_postprocessor.cc new file mode 100644 index 0000000000..e018456c4a --- /dev/null +++ b/deal.II/deal.II/source/numerics/data_postprocessor.cc @@ -0,0 +1,53 @@ +//--------------------------------------------------------------------------- +// $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 + +DEAL_II_NAMESPACE_OPEN + + +template +DataPostprocessor::~DataPostprocessor() +{} + + + +template +void +DataPostprocessor::compute_derived_quantities (std::vector > &computed_quantities, + const std::vector &/*uh*/, + const std::vector > &/*duh*/, + const std::vector > &/*dduh*/, + const std::vector > &/*normals*/) const +{ + computed_quantities.clear(); + AssertThrow(false,ExcPureFunctionCalled()); +} + + + +template +void +DataPostprocessor::compute_derived_quantities (std::vector > &computed_quantities, + const std::vector > &/*uh*/, + const std::vector > > &/*duh*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/) const +{ + computed_quantities.clear(); + AssertThrow(false,ExcPureFunctionCalled()); +} + +// explicit instantiation +template class DataPostprocessor; + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5