From d01d4cefffa9ec6bea169661ddb1042afad2f800 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 13 May 2017 15:56:26 -0400 Subject: [PATCH] Split DataOut_DoFData instantiations. Before this commit GCC needed 2.3 GB, 65s in release mode and 3.6 GB, 79s in debug mode to compile this file. After splitting the file we need 43s/32s and 1.6 GB/1.3 GB in release mode and 52s/34s and 2.1 GB/2.0 GB in debug mode. This is substantial because this was (previously) the most expensive file to compile. This also adds a few missing instantiations to the nonzero codimension case. --- .../numerics/data_out_dof_data.templates.h | 1281 +++++++++++++++++ source/numerics/CMakeLists.txt | 2 + source/numerics/data_out_dof_data.cc | 1262 +--------------- source/numerics/data_out_dof_data.inst.in | 94 +- source/numerics/data_out_dof_data_codim.cc | 22 + .../numerics/data_out_dof_data_codim.inst.in | 106 ++ 6 files changed, 1424 insertions(+), 1343 deletions(-) create mode 100644 include/deal.II/numerics/data_out_dof_data.templates.h create mode 100644 source/numerics/data_out_dof_data_codim.cc create mode 100644 source/numerics/data_out_dof_data_codim.inst.in diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h new file mode 100644 index 0000000000..18ff0a6624 --- /dev/null +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -0,0 +1,1281 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1999 - 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__data_out_dof_data_templates_h +#define dealii__data_out_dof_data_templates_h + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + namespace DataOut + { + template + ParallelDataBase:: + ParallelDataBase (const unsigned int n_datasets, + const unsigned int n_subdivisions, + const std::vector &n_postprocessor_outputs, + const Mapping &mapping, + const std::vector > > &finite_elements, + const UpdateFlags update_flags, + const bool use_face_values) + : + n_datasets (n_datasets), + n_subdivisions (n_subdivisions), + postprocessed_values (n_postprocessor_outputs.size()), + mapping_collection (mapping), + finite_elements (finite_elements), + update_flags (update_flags) + { + unsigned int n_q_points = 0; + if (use_face_values == false) + { + dealii::hp::QCollection + quadrature(QIterated(QTrapez<1>(), n_subdivisions)); + n_q_points = quadrature[0].size(); + x_fe_values.resize(this->finite_elements.size()); + for (unsigned int i=0; ifinite_elements.size(); ++i) + { + // check if there is a finite element that is equal to the + // present one, then we can re-use the FEValues object + for (unsigned int j=0; jfinite_elements[i].get() == + this->finite_elements[j].get()) + { + x_fe_values[i] = x_fe_values[j]; + break; + } + if (x_fe_values[i].get() == nullptr) + x_fe_values[i].reset(new dealii::hp::FEValues + (this->mapping_collection, + *this->finite_elements[i], + quadrature, + this->update_flags)); + } + } + else + { + dealii::hp::QCollection + quadrature(QIterated(QTrapez<1>(), n_subdivisions)); + n_q_points = quadrature[0].size(); + x_fe_face_values.resize(this->finite_elements.size()); + for (unsigned int i=0; ifinite_elements.size(); ++i) + { + // check if there is a finite element that is equal to the + // present one, then we can re-use the FEValues object + for (unsigned int j=0; jfinite_elements[i].get() == + this->finite_elements[j].get()) + { + x_fe_face_values[i] = x_fe_face_values[j]; + break; + } + if (x_fe_face_values[i].get() == nullptr) + x_fe_face_values[i].reset(new dealii::hp::FEFaceValues + (this->mapping_collection, + *this->finite_elements[i], + quadrature, + this->update_flags)); + } + } + + patch_values_scalar.solution_values.resize (n_q_points); + patch_values_scalar.solution_gradients.resize (n_q_points); + patch_values_scalar.solution_hessians.resize (n_q_points); + patch_values_system.solution_values.resize (n_q_points); + patch_values_system.solution_gradients.resize (n_q_points); + patch_values_system.solution_hessians.resize (n_q_points); + + for (unsigned int dataset=0; dataset(n_postprocessor_outputs[dataset])); + } + + + + + + // implement copy constructor to create a thread's own version of + // x_fe_values + template + ParallelDataBase:: + ParallelDataBase (const ParallelDataBase &data) + : + n_datasets (data.n_datasets), + n_subdivisions (data.n_subdivisions), + patch_values_scalar (data.patch_values_scalar), + patch_values_system (data.patch_values_system), + postprocessed_values (data.postprocessed_values), + mapping_collection (data.mapping_collection), + finite_elements (data.finite_elements), + update_flags (data.update_flags) + { + if (data.x_fe_values.empty() == false) + { + Assert(data.x_fe_face_values.empty() == true, ExcInternalError()); + dealii::hp::QCollection + quadrature(QIterated(QTrapez<1>(), n_subdivisions)); + x_fe_values.resize(this->finite_elements.size()); + for (unsigned int i=0; ifinite_elements.size(); ++i) + { + // check if there is a finite element that is equal to the + // present one, then we can re-use the FEValues object + for (unsigned int j=0; jfinite_elements[i].get() == + this->finite_elements[j].get()) + { + x_fe_values[i] = x_fe_values[j]; + break; + } + if (x_fe_values[i].get() == nullptr) + x_fe_values[i].reset(new dealii::hp::FEValues + (this->mapping_collection, + *this->finite_elements[i], + quadrature, + this->update_flags)); + } + } + else + { + dealii::hp::QCollection + quadrature(QIterated(QTrapez<1>(), n_subdivisions)); + x_fe_face_values.resize(this->finite_elements.size()); + for (unsigned int i=0; ifinite_elements.size(); ++i) + { + // check if there is a finite element that is equal to the + // present one, then we can re-use the FEValues object + for (unsigned int j=0; jfinite_elements[i].get() == + this->finite_elements[j].get()) + { + x_fe_face_values[i] = x_fe_face_values[j]; + break; + } + if (x_fe_face_values[i].get() == nullptr) + x_fe_face_values[i].reset(new dealii::hp::FEFaceValues + (this->mapping_collection, + *this->finite_elements[i], + quadrature, + this->update_flags)); + } + } + } + + + + template + template + void + ParallelDataBase:: + reinit_all_fe_values(std::vector > > &dof_data, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face) + { + for (unsigned int dataset=0; datasetget_triangulation(), + cell->level(), + cell->index(), + dof_data[dataset]->dof_handler); + if (x_fe_values.empty()) + { + AssertIndexRange(face, + GeometryInfo::faces_per_cell); + x_fe_face_values[dataset]->reinit(dh_cell, face); + } + else + x_fe_values[dataset]->reinit (dh_cell); + } + } + if (dof_data.empty()) + { + if (x_fe_values.empty()) + { + AssertIndexRange(face, + GeometryInfo::faces_per_cell); + x_fe_face_values[0]->reinit(cell, face); + } + else + x_fe_values[0]->reinit (cell); + } + } + + + + template + const FEValuesBase & + ParallelDataBase:: + get_present_fe_values(const unsigned int dataset) const + { + AssertIndexRange(dataset, finite_elements.size()); + if (x_fe_values.empty()) + return x_fe_face_values[dataset]->get_present_fe_values(); + else + return x_fe_values[dataset]->get_present_fe_values(); + } + + + + template + void + ParallelDataBase:: + resize_system_vectors(const unsigned int n_components) + { + Assert(patch_values_system.solution_values.size() > 0, ExcInternalError()); + AssertDimension(patch_values_system.solution_values.size(), + patch_values_system.solution_gradients.size()); + AssertDimension(patch_values_system.solution_values.size(), + patch_values_system.solution_hessians.size()); + if (patch_values_system.solution_values[0].size() == n_components) + return; + for (unsigned int k=0; k + void + append_patch_to_list (const DataOutBase::Patch &patch, + std::vector > &patches) + { + patches.push_back (patch); + patches.back().patch_index = patches.size()-1; + } + } +} + +namespace internal +{ + namespace DataOut + { + template + DataEntryBase::DataEntryBase + (const DoFHandlerType *dofs, + const std::vector &names_in, + const std::vector &data_component_interpretation) + : + dof_handler (dofs, typeid(dealii::DataOut_DoFData).name()), + names(names_in), + data_component_interpretation (data_component_interpretation), + postprocessor(nullptr, typeid(*this).name()), + n_output_variables(names.size()) + { + Assert (names.size() == data_component_interpretation.size(), + ExcDimensionMismatch(data_component_interpretation.size(), + names.size())); + + // check that the names use only allowed characters + for (unsigned int i=0; i()") == std::string::npos, + Exceptions::DataOut::ExcInvalidCharacter (names[i], + names[i].find_first_not_of("abcdefghijklmnopqrstuvwxyz" + "ABCDEFGHIJKLMNOPQRSTUVWXYZ" + "0123456789_<>()"))); + } + + + + template + DataEntryBase::DataEntryBase + (const DoFHandlerType *dofs, + const DataPostprocessor *data_postprocessor) + : + dof_handler (dofs, typeid(dealii::DataOut_DoFData).name()), + names(data_postprocessor->get_names()), + data_component_interpretation (data_postprocessor->get_data_component_interpretation()), + postprocessor(data_postprocessor, typeid(*this).name()), + n_output_variables(names.size()) + { + Assert (data_postprocessor->get_names().size() + == + data_postprocessor->get_data_component_interpretation().size(), + ExcDimensionMismatch (data_postprocessor->get_names().size(), + data_postprocessor->get_data_component_interpretation().size())); + + // check that the names use only allowed characters + for (unsigned int i=0; i()") == std::string::npos, + Exceptions::DataOut::ExcInvalidCharacter (names[i], + names[i].find_first_not_of("abcdefghijklmnopqrstuvwxyz" + "ABCDEFGHIJKLMNOPQRSTUVWXYZ" + "0123456789_<>()"))); + } + + + + template + DataEntryBase::~DataEntryBase () + {} + + + + /** + * Class that stores a pointer to a vector of type equal to the template + * argument, and provides the functions to extract data from it. + * + * @author Wolfgang Bangerth, 2004 + */ + template + class DataEntry : public DataEntryBase + { + public: + /** + * Constructor. Give a list of names for the individual components of + * the vector and their interpretation as scalar or vector data. This + * constructor assumes that no postprocessor is going to be used. + */ + DataEntry + (const DoFHandlerType *dofs, + const VectorType *data, + const std::vector &names, + const std::vector &data_component_interpretation); + + /** + * Constructor when a data postprocessor is going to be used. In that + * case, the names and vector declarations are going to be acquired from + * the postprocessor. + */ + DataEntry (const DoFHandlerType *dofs, + const VectorType *data, + const DataPostprocessor *data_postprocessor); + + /** + * Assuming that the stored vector is a cell vector, extract the given + * element from it. + */ + virtual + double + get_cell_data_value (const unsigned int cell_number) const; + + /** + * Given a FEValuesBase object, extract the values on the present cell + * from the vector we actually store. + */ + virtual + void + get_function_values + (const FEValuesBase &fe_patch_values, + std::vector &patch_values) const; + + /** + * Given a FEValuesBase object, extract the values 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_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_hessians + (const FEValuesBase &fe_patch_values, + std::vector > &patch_hessians) 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_hessians + (const FEValuesBase &fe_patch_values, + std::vector > > &patch_hessians_system) const; + + /** + * Clear all references to the vectors. + */ + virtual void clear (); + + /** + * Determine an estimate for the memory consumption (in bytes) of this + * object. + */ + virtual std::size_t memory_consumption () const; + + private: + /** + * Pointer to the data vector. Note that ownership of the vector pointed + * to remains with the caller of this class. + */ + const VectorType *vector; + }; + + + + template + DataEntry:: + DataEntry (const DoFHandlerType *dofs, + const VectorType *data, + const std::vector &names, + const std::vector &data_component_interpretation) + : + DataEntryBase (dofs, names, data_component_interpretation), + vector (data) + {} + + + + template + DataEntry:: + DataEntry (const DoFHandlerType *dofs, + const VectorType *data, + const DataPostprocessor *data_postprocessor) + : + DataEntryBase (dofs, data_postprocessor), + vector (data) + {} + + + namespace + { + template + double + get_vector_element (const VectorType &vector, + const unsigned int cell_number) + { + return internal::ElementAccess::get(vector,cell_number); + } + + + inline double + get_vector_element (const IndexSet &is, + const unsigned int cell_number) + { + return (is.is_element(cell_number) ? 1 : 0); + } + } + + + + template + double + DataEntry:: + get_cell_data_value (const unsigned int cell_number) const + { + return get_vector_element(*vector, cell_number); + } + + + + template + void + DataEntry::get_function_values + (const FEValuesBase &fe_patch_values, + std::vector > &patch_values_system) const + { + // FIXME: FEValuesBase gives us data in types that match that of + // the solution vector. but this function needs to pass it back + // up as 'double' vectors. this requires the use of a temporary + // variable here if the data we get is not a 'double' vector. + // (of course, in reality, this also means that we may lose + // information to begin with.) + // + // the correct thing would be to also use the correct data type + // upstream somewhere, but this is complicated because we hide + // the actual data type from upstream. rather, we should at + // least make sure we can deal with complex numbers + if (typeid(typename VectorType::value_type) == typeid(double)) + { + fe_patch_values.get_function_values (*vector, + // reinterpret output argument type; because of + // the 'if' statement above, this is the + // identity cast whenever the code is + // executed, but the cast is necessary + // to allow compilation even if we don't get here + reinterpret_cast >&> + (patch_values_system)); + } + else + { + std::vector > tmp(patch_values_system.size()); + for (unsigned int i = 0; i < patch_values_system.size(); i++) + tmp[i].reinit(patch_values_system[i]); + + fe_patch_values.get_function_values (*vector, tmp); + + for (unsigned int i = 0; i < patch_values_system.size(); i++) + patch_values_system[i] = tmp[i]; + } + } + + + + template + void + DataEntry::get_function_values + (const FEValuesBase &fe_patch_values, + std::vector &patch_values) const + { + // FIXME: FEValuesBase gives us data in types that match that of + // the solution vector. but this function needs to pass it back + // up as 'double' vectors. this requires the use of a temporary + // variable here if the data we get is not a 'double' vector. + // (of course, in reality, this also means that we may lose + // information to begin with.) + // + // the correct thing would be to also use the correct data type + // upstream somewhere, but this is complicated because we hide + // the actual data type from upstream. rather, we should at + // least make sure we can deal with complex numbers + if (typeid(typename VectorType::value_type) == typeid(double)) + { + fe_patch_values.get_function_values (*vector, + // reinterpret output argument type; because of + // the 'if' statement above, this is the + // identity cast whenever the code is + // executed, but the cast is necessary + // to allow compilation even if we don't get here + reinterpret_cast&> + (patch_values)); + } + else + { + std::vector tmp (patch_values.size()); + + fe_patch_values.get_function_values (*vector, tmp); + + for (unsigned int i = 0; i < tmp.size(); i++) + patch_values[i] = tmp[i]; + } + } + + + + template + void + DataEntry::get_function_gradients + (const FEValuesBase &fe_patch_values, + std::vector > > &patch_gradients_system) const + { + // FIXME: FEValuesBase gives us data in types that match that of + // the solution vector. but this function needs to pass it back + // up as 'double' vectors. this requires the use of a temporary + // variable here if the data we get is not a 'double' vector. + // (of course, in reality, this also means that we may lose + // information to begin with.) + // + // the correct thing would be to also use the correct data type + // upstream somewhere, but this is complicated because we hide + // the actual data type from upstream. rather, we should at + // least make sure we can deal with complex numbers + if (typeid(typename VectorType::value_type) == typeid(double)) + { + fe_patch_values.get_function_gradients (*vector, + // reinterpret output argument type; because of + // the 'if' statement above, this is the + // identity cast whenever the code is + // executed, but the cast is necessary + // to allow compilation even if we don't get here + reinterpret_cast > >&> + (patch_gradients_system)); + } + else + { + std::vector > > + tmp(patch_gradients_system.size()); + for (unsigned int i = 0; i < tmp.size(); i++) + tmp[i].resize(patch_gradients_system[i].size()); + + fe_patch_values.get_function_gradients (*vector, tmp); + + for (unsigned int i = 0; i < tmp.size(); i++) + for (unsigned int j = 0; j < tmp[i].size(); j++) + patch_gradients_system[i][j] = tmp[i][j]; + } + } + + + + template + void + DataEntry::get_function_gradients + (const FEValuesBase &fe_patch_values, + std::vector > &patch_gradients) const + { + // FIXME: FEValuesBase gives us data in types that match that of + // the solution vector. but this function needs to pass it back + // up as 'double' vectors. this requires the use of a temporary + // variable here if the data we get is not a 'double' vector. + // (of course, in reality, this also means that we may lose + // information to begin with.) + // + // the correct thing would be to also use the correct data type + // upstream somewhere, but this is complicated because we hide + // the actual data type from upstream. rather, we should at + // least make sure we can deal with complex numbers + if (typeid(typename VectorType::value_type) == typeid(double)) + { + fe_patch_values.get_function_gradients (*vector, + // reinterpret output argument type; because of + // the 'if' statement above, this is the + // identity cast whenever the code is + // executed, but the cast is necessary + // to allow compilation even if we don't get here + reinterpret_cast >&> + (patch_gradients)); + } + else + { + std::vector > tmp; + tmp.resize(patch_gradients.size()); + + fe_patch_values.get_function_gradients (*vector, tmp); + + for (unsigned int i = 0; i < tmp.size(); i++) + patch_gradients[i] = tmp[i]; + } + } + + + + template + void + DataEntry::get_function_hessians + (const FEValuesBase &fe_patch_values, + std::vector > > &patch_hessians_system) const + { + // FIXME: FEValuesBase gives us data in types that match that of + // the solution vector. but this function needs to pass it back + // up as 'double' vectors. this requires the use of a temporary + // variable here if the data we get is not a 'double' vector. + // (of course, in reality, this also means that we may lose + // information to begin with.) + // + // the correct thing would be to also use the correct data type + // upstream somewhere, but this is complicated because we hide + // the actual data type from upstream. rather, we should at + // least make sure we can deal with complex numbers + if (typeid(typename VectorType::value_type) == typeid(double)) + { + fe_patch_values.get_function_hessians (*vector, + // reinterpret output argument type; because of + // the 'if' statement above, this is the + // identity cast whenever the code is + // executed, but the cast is necessary + // to allow compilation even if we don't get here + reinterpret_cast > >&> + (patch_hessians_system)); + } + else + { + std::vector > > + tmp(patch_hessians_system.size()); + for (unsigned int i = 0; i < tmp.size(); i++) + tmp[i].resize(patch_hessians_system[i].size()); + + fe_patch_values.get_function_hessians (*vector, tmp); + + for (unsigned int i = 0; i < tmp.size(); i++) + for (unsigned int j = 0; j < tmp[i].size(); j++) + patch_hessians_system[i][j] = tmp[i][j]; + } + } + + + + template + void + DataEntry::get_function_hessians + (const FEValuesBase &fe_patch_values, + std::vector > &patch_hessians) const + { + // FIXME: FEValuesBase gives us data in types that match that of + // the solution vector. but this function needs to pass it back + // up as 'double' vectors. this requires the use of a temporary + // variable here if the data we get is not a 'double' vector. + // (of course, in reality, this also means that we may lose + // information to begin with.) + // + // the correct thing would be to also use the correct data type + // upstream somewhere, but this is complicated because we hide + // the actual data type from upstream. rather, we should at + // least make sure we can deal with complex numbers + if (typeid(typename VectorType::value_type) == typeid(double)) + { + fe_patch_values.get_function_hessians (*vector, + // reinterpret output argument type; because of + // the 'if' statement above, this is the + // identity cast whenever the code is + // executed, but the cast is necessary + // to allow compilation even if we don't get here + reinterpret_cast >&> + (patch_hessians)); + } + else + { + std::vector > + tmp(patch_hessians.size()); + + fe_patch_values.get_function_hessians (*vector, tmp); + + for (unsigned int i = 0; i < tmp.size(); i++) + patch_hessians[i] = tmp[i]; + } + } + + + + template + std::size_t + DataEntry::memory_consumption () const + { + return (sizeof (vector) + + MemoryConsumption::memory_consumption (this->names)); + } + + + + template + void + DataEntry::clear () + { + vector = nullptr; + this->dof_handler = nullptr; + } + } +} + + + +template +DataOut_DoFData::DataOut_DoFData () + : + triangulation(nullptr,typeid(*this).name()), + dofs(nullptr,typeid(*this).name()) +{} + + + +template +DataOut_DoFData::~DataOut_DoFData () +{ + clear (); +} + + + +template +void +DataOut_DoFData:: +attach_dof_handler (const DoFHandlerType &d) +{ + Assert (dof_data.size() == 0, + Exceptions::DataOut::ExcOldDataStillPresent()); + Assert (cell_data.size() == 0, + Exceptions::DataOut::ExcOldDataStillPresent()); + + triangulation = SmartPointer > + (&d.get_triangulation(), typeid(*this).name()); + dofs = SmartPointer(&d, typeid(*this).name()); +} + + + +template +void +DataOut_DoFData:: +attach_triangulation (const Triangulation &tria) +{ + Assert (dof_data.size() == 0, + Exceptions::DataOut::ExcOldDataStillPresent()); + Assert (cell_data.size() == 0, + Exceptions::DataOut::ExcOldDataStillPresent()); + + triangulation = SmartPointer > + (&tria, typeid(*this).name()); +} + + + +template +template +void +DataOut_DoFData:: +add_data_vector (const DoFHandlerType &dof_handler, + const VectorType &vec, + const DataPostprocessor &data_postprocessor) +{ + // this is a specialized version of the other function where we have a + // postprocessor. if we do, we know that we have type_dof_data, which makes + // things a bit simpler, we also don't need to deal with some of the other + // stuff and use a different constructor of DataEntry + if (triangulation != nullptr) + { + Assert (&dof_handler.get_triangulation() == triangulation, + ExcMessage("The triangulation attached to the DoFHandler does not " + "match with the one set previously")); + } + else + { + triangulation = SmartPointer> + (&dof_handler.get_triangulation(), typeid(*this).name()); + } + + Assert (vec.size() == dof_handler.n_dofs(), + Exceptions::DataOut::ExcInvalidVectorSize (vec.size(), + dof_handler.n_dofs(), + dof_handler.get_triangulation().n_active_cells())); + + + auto new_entry = std_cxx14::make_unique> + (&dof_handler, &vec, &data_postprocessor); + dof_data.emplace_back (std::move(new_entry)); +} + + + +template +template +void +DataOut_DoFData:: +add_data_vector_internal +(const DoFHandlerType *dof_handler, + const VectorType &data_vector, + const std::vector &names, + const DataVectorType type, + const std::vector &data_component_interpretation_, + const bool deduce_output_names) +{ + // Check available mesh information: + if (triangulation == nullptr) + { + Assert(dof_handler != nullptr, ExcInternalError()); + triangulation = SmartPointer> + (&dof_handler->get_triangulation(), typeid(*this).name()); + } + + if (dof_handler != nullptr) + { + Assert (&dof_handler->get_triangulation() == triangulation, + ExcMessage("The triangulation attached to the DoFHandler does not " + "match with the one set previously")); + } + + // Figure out the data type: + DataVectorType actual_type = type; + if (type == type_automatic) + { + Assert((dof_handler == nullptr) || (triangulation->n_active_cells() != dof_handler->n_dofs()), + ExcMessage("Unable to determine the type of vector automatically because the number of DoFs " + "is equal to the number of cells. Please specify DataVectorType.")); + + if (data_vector.size() == triangulation->n_active_cells()) + actual_type = type_cell_data; + else + actual_type = type_dof_data; + } + Assert(actual_type == type_cell_data || actual_type == type_dof_data, + ExcInternalError()); + + // If necessary, append '_1', '_2', etc. to component names: + std::vector deduced_names; + if (deduce_output_names && actual_type == type_dof_data) + { + Assert(names.size() == 1, ExcInternalError()); + Assert(dof_handler != nullptr, ExcInternalError()); + Assert(dof_handler->n_dofs() > 0, + ExcMessage("The DoF handler attached to the current output vector " + "does not have any degrees of freedom, so it is not " + "possible to output DoF data in this context.")); + const std::string name = names[0]; + const unsigned int n_components = dof_handler->get_fe().n_components(); + deduced_names.resize (n_components); + if (n_components > 1) + { + for (unsigned int i=0; in_active_cells(), + ExcDimensionMismatch (data_vector.size(), + triangulation->n_active_cells())); + Assert (deduced_names.size() == 1, + Exceptions::DataOut::ExcInvalidNumberOfNames (deduced_names.size(), 1)); + break; + case type_dof_data: + Assert (dof_handler != nullptr, + Exceptions::DataOut::ExcNoDoFHandlerSelected ()); + Assert (data_vector.size() == dof_handler->n_dofs(), + Exceptions::DataOut::ExcInvalidVectorSize (data_vector.size(), + dof_handler->n_dofs(), + triangulation->n_active_cells())); + Assert (deduced_names.size() == dof_handler->get_fe().n_components(), + Exceptions::DataOut::ExcInvalidNumberOfNames (deduced_names.size(), + dof_handler->get_fe().n_components())); + break; + default: + Assert (false, ExcInternalError()); + } + + const auto &data_component_interpretation + = (data_component_interpretation_.size() != 0 + ? + data_component_interpretation_ + : + std::vector + (deduced_names.size(), DataComponentInterpretation::component_is_scalar)); + + // finally, add the data vector: + auto new_entry = std_cxx14::make_unique> + (dof_handler, &data_vector, deduced_names, data_component_interpretation); + + if (actual_type == type_dof_data) + dof_data.emplace_back (std::move(new_entry)); + else + cell_data.emplace_back (std::move(new_entry)); +} + + + +template +void DataOut_DoFData::clear_data_vectors () +{ + dof_data.erase (dof_data.begin(), dof_data.end()); + cell_data.erase (cell_data.begin(), cell_data.end()); + + // delete patches + std::vector dummy; + patches.swap (dummy); +} + + + +template +void +DataOut_DoFData:: +clear_input_data_references () +{ + for (unsigned int i=0; iclear (); + + for (unsigned int i=0; iclear (); + + if (dofs != nullptr) + dofs = nullptr; +} + + + +template +void +DataOut_DoFData::clear () +{ + dof_data.erase (dof_data.begin(), dof_data.end()); + cell_data.erase (cell_data.begin(), cell_data.end()); + + if (dofs != nullptr) + dofs = nullptr; + + // delete patches + std::vector dummy; + patches.swap (dummy); +} + + + +template +std::vector +DataOut_DoFData:: +get_dataset_names () const +{ + std::vector names; + // collect the names of dof + // and cell data + typedef + typename std::vector > >::const_iterator + data_iterator; + + for (data_iterator d=dof_data.begin(); + d!=dof_data.end(); ++d) + for (unsigned int i=0; i<(*d)->names.size(); ++i) + names.push_back ((*d)->names[i]); + for (data_iterator d=cell_data.begin(); d!=cell_data.end(); ++d) + { + Assert ((*d)->names.size() == 1, ExcInternalError()); + names.push_back ((*d)->names[0]); + } + + return names; +} + + + +template +std::vector > +DataOut_DoFData::get_vector_data_ranges () const +{ + std::vector > + ranges; + + // collect the ranges of dof + // and cell data + typedef + typename std::vector > >::const_iterator + data_iterator; + + unsigned int output_component = 0; + for (data_iterator d=dof_data.begin(); + d!=dof_data.end(); ++d) + for (unsigned int i=0; i<(*d)->n_output_variables; + ++i, ++output_component) + // see what kind of data we have + // here. note that for the purpose of + // the current function all we care + // about is vector data + if ((*d)->data_component_interpretation[i] == + DataComponentInterpretation::component_is_part_of_vector) + { + // ensure that there is a + // continuous number of next + // space_dim components that all + // deal with vectors + Assert (i+patch_space_dim <= + (*d)->n_output_variables, + Exceptions::DataOut::ExcInvalidVectorDeclaration (i, + (*d)->names[i])); + for (unsigned int dd=1; dddata_component_interpretation[i+dd] + == + DataComponentInterpretation::component_is_part_of_vector, + Exceptions::DataOut::ExcInvalidVectorDeclaration (i, + (*d)->names[i])); + + // all seems alright, so figure out + // whether there is a common name + // to these components. if not, + // leave the name empty and let the + // output format writer decide what + // to do here + std::string name = (*d)->names[i]; + for (unsigned int dd=1; ddnames[i+dd]) + { + name = ""; + break; + } + + // finally add a corresponding + // range + std::tuple + range (output_component, + output_component+patch_space_dim-1, + name); + + ranges.push_back (range); + + // increase the 'component' counter + // by the appropriate amount, same + // for 'i', since we have already + // dealt with all these components + output_component += patch_space_dim-1; + i += patch_space_dim-1; + } + + // note that we do not have to traverse the + // list of cell data here because cell data + // is one value per (logical) cell and + // therefore cannot be a vector + + // as a final check, the 'component' + // counter should be at the total number of + // components added up now +#ifdef DEBUG + unsigned int n_output_components = 0; + for (data_iterator d=dof_data.begin(); + d!=dof_data.end(); ++d) + n_output_components += (*d)->n_output_variables; + Assert (output_component == n_output_components, + ExcInternalError()); +#endif + + return ranges; +} + + + +template +const std::vector< dealii::DataOutBase::Patch > & +DataOut_DoFData::get_patches () const +{ + return patches; +} + + + +template +std::vector > > + DataOut_DoFData::get_finite_elements() const +{ + const unsigned int dhdim = DoFHandlerType::dimension; + const unsigned int dhspacedim = DoFHandlerType::space_dimension; + std::vector > > + finite_elements(this->dof_data.size()); + for (unsigned int i=0; idof_data.size(); ++i) + { + Assert (dof_data[i]->dof_handler != nullptr, + Exceptions::DataOut::ExcNoDoFHandlerSelected ()); + + // avoid creating too many finite elements and doing a lot of work on + // initializing FEValues downstream: if two DoFHandlers are the same + // (checked by pointer comparison), we can re-use the shared_ptr object + // for the second one. We cannot check for finite element equalities + // because we need different FEValues objects for different dof + // handlers. + bool duplicate = false; + for (unsigned int j=0; jdof_handler == dof_data[j]->dof_handler) + { + finite_elements[i] = finite_elements[j]; + duplicate = true; + } + if (duplicate == false) + finite_elements[i].reset(new dealii::hp::FECollection + (this->dof_data[i]->dof_handler->get_fe())); + } + if (this->dof_data.empty()) + { + finite_elements.resize(1); + finite_elements[0].reset(new dealii::hp::FECollection + (FE_DGQ(0))); + } + return finite_elements; +} + + + +template +std::size_t +DataOut_DoFData::memory_consumption () const +{ + return (DataOutInterface::memory_consumption () + + MemoryConsumption::memory_consumption (dofs) + + MemoryConsumption::memory_consumption (patches)); +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/numerics/CMakeLists.txt b/source/numerics/CMakeLists.txt index 9f1fb1f340..85bbf2fabc 100644 --- a/source/numerics/CMakeLists.txt +++ b/source/numerics/CMakeLists.txt @@ -18,6 +18,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src data_out.cc data_out_dof_data.cc + data_out_dof_data_codim.cc data_out_faces.cc data_out_rotation.cc data_out_stack.cc @@ -59,6 +60,7 @@ SET(_src SET(_inst data_out_dof_data.inst.in + data_out_dof_data_codim.inst.in data_out_faces.inst.in data_out.inst.in data_out_rotation.inst.in diff --git a/source/numerics/data_out_dof_data.cc b/source/numerics/data_out_dof_data.cc index 00d89adee2..5ec70d1eec 100644 --- a/source/numerics/data_out_dof_data.cc +++ b/source/numerics/data_out_dof_data.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 1999 - 2016 by the deal.II authors +// Copyright (C) 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -13,1268 +13,10 @@ // // --------------------------------------------------------------------- -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include +#include DEAL_II_NAMESPACE_OPEN - -namespace internal -{ - namespace DataOut - { - template - ParallelDataBase:: - ParallelDataBase (const unsigned int n_datasets, - const unsigned int n_subdivisions, - const std::vector &n_postprocessor_outputs, - const Mapping &mapping, - const std::vector > > &finite_elements, - const UpdateFlags update_flags, - const bool use_face_values) - : - n_datasets (n_datasets), - n_subdivisions (n_subdivisions), - postprocessed_values (n_postprocessor_outputs.size()), - mapping_collection (mapping), - finite_elements (finite_elements), - update_flags (update_flags) - { - unsigned int n_q_points = 0; - if (use_face_values == false) - { - dealii::hp::QCollection - quadrature(QIterated(QTrapez<1>(), n_subdivisions)); - n_q_points = quadrature[0].size(); - x_fe_values.resize(this->finite_elements.size()); - for (unsigned int i=0; ifinite_elements.size(); ++i) - { - // check if there is a finite element that is equal to the - // present one, then we can re-use the FEValues object - for (unsigned int j=0; jfinite_elements[i].get() == - this->finite_elements[j].get()) - { - x_fe_values[i] = x_fe_values[j]; - break; - } - if (x_fe_values[i].get() == nullptr) - x_fe_values[i].reset(new dealii::hp::FEValues - (this->mapping_collection, - *this->finite_elements[i], - quadrature, - this->update_flags)); - } - } - else - { - dealii::hp::QCollection - quadrature(QIterated(QTrapez<1>(), n_subdivisions)); - n_q_points = quadrature[0].size(); - x_fe_face_values.resize(this->finite_elements.size()); - for (unsigned int i=0; ifinite_elements.size(); ++i) - { - // check if there is a finite element that is equal to the - // present one, then we can re-use the FEValues object - for (unsigned int j=0; jfinite_elements[i].get() == - this->finite_elements[j].get()) - { - x_fe_face_values[i] = x_fe_face_values[j]; - break; - } - if (x_fe_face_values[i].get() == nullptr) - x_fe_face_values[i].reset(new dealii::hp::FEFaceValues - (this->mapping_collection, - *this->finite_elements[i], - quadrature, - this->update_flags)); - } - } - - patch_values_scalar.solution_values.resize (n_q_points); - patch_values_scalar.solution_gradients.resize (n_q_points); - patch_values_scalar.solution_hessians.resize (n_q_points); - patch_values_system.solution_values.resize (n_q_points); - patch_values_system.solution_gradients.resize (n_q_points); - patch_values_system.solution_hessians.resize (n_q_points); - - for (unsigned int dataset=0; dataset(n_postprocessor_outputs[dataset])); - } - - - - - - // implement copy constructor to create a thread's own version of - // x_fe_values - template - ParallelDataBase:: - ParallelDataBase (const ParallelDataBase &data) - : - n_datasets (data.n_datasets), - n_subdivisions (data.n_subdivisions), - patch_values_scalar (data.patch_values_scalar), - patch_values_system (data.patch_values_system), - postprocessed_values (data.postprocessed_values), - mapping_collection (data.mapping_collection), - finite_elements (data.finite_elements), - update_flags (data.update_flags) - { - if (data.x_fe_values.empty() == false) - { - Assert(data.x_fe_face_values.empty() == true, ExcInternalError()); - dealii::hp::QCollection - quadrature(QIterated(QTrapez<1>(), n_subdivisions)); - x_fe_values.resize(this->finite_elements.size()); - for (unsigned int i=0; ifinite_elements.size(); ++i) - { - // check if there is a finite element that is equal to the - // present one, then we can re-use the FEValues object - for (unsigned int j=0; jfinite_elements[i].get() == - this->finite_elements[j].get()) - { - x_fe_values[i] = x_fe_values[j]; - break; - } - if (x_fe_values[i].get() == nullptr) - x_fe_values[i].reset(new dealii::hp::FEValues - (this->mapping_collection, - *this->finite_elements[i], - quadrature, - this->update_flags)); - } - } - else - { - dealii::hp::QCollection - quadrature(QIterated(QTrapez<1>(), n_subdivisions)); - x_fe_face_values.resize(this->finite_elements.size()); - for (unsigned int i=0; ifinite_elements.size(); ++i) - { - // check if there is a finite element that is equal to the - // present one, then we can re-use the FEValues object - for (unsigned int j=0; jfinite_elements[i].get() == - this->finite_elements[j].get()) - { - x_fe_face_values[i] = x_fe_face_values[j]; - break; - } - if (x_fe_face_values[i].get() == nullptr) - x_fe_face_values[i].reset(new dealii::hp::FEFaceValues - (this->mapping_collection, - *this->finite_elements[i], - quadrature, - this->update_flags)); - } - } - } - - - - template - template - void - ParallelDataBase:: - reinit_all_fe_values(std::vector > > &dof_data, - const typename dealii::Triangulation::cell_iterator &cell, - const unsigned int face) - { - for (unsigned int dataset=0; datasetget_triangulation(), - cell->level(), - cell->index(), - dof_data[dataset]->dof_handler); - if (x_fe_values.empty()) - { - AssertIndexRange(face, - GeometryInfo::faces_per_cell); - x_fe_face_values[dataset]->reinit(dh_cell, face); - } - else - x_fe_values[dataset]->reinit (dh_cell); - } - } - if (dof_data.empty()) - { - if (x_fe_values.empty()) - { - AssertIndexRange(face, - GeometryInfo::faces_per_cell); - x_fe_face_values[0]->reinit(cell, face); - } - else - x_fe_values[0]->reinit (cell); - } - } - - - - template - const FEValuesBase & - ParallelDataBase:: - get_present_fe_values(const unsigned int dataset) const - { - AssertIndexRange(dataset, finite_elements.size()); - if (x_fe_values.empty()) - return x_fe_face_values[dataset]->get_present_fe_values(); - else - return x_fe_values[dataset]->get_present_fe_values(); - } - - - - template - void - ParallelDataBase:: - resize_system_vectors(const unsigned int n_components) - { - Assert(patch_values_system.solution_values.size() > 0, ExcInternalError()); - AssertDimension(patch_values_system.solution_values.size(), - patch_values_system.solution_gradients.size()); - AssertDimension(patch_values_system.solution_values.size(), - patch_values_system.solution_hessians.size()); - if (patch_values_system.solution_values[0].size() == n_components) - return; - for (unsigned int k=0; k - void - append_patch_to_list (const DataOutBase::Patch &patch, - std::vector > &patches) - { - patches.push_back (patch); - patches.back().patch_index = patches.size()-1; - } - } -} - -namespace internal -{ - namespace DataOut - { - template - DataEntryBase::DataEntryBase - (const DoFHandlerType *dofs, - const std::vector &names_in, - const std::vector &data_component_interpretation) - : - dof_handler (dofs, typeid(dealii::DataOut_DoFData).name()), - names(names_in), - data_component_interpretation (data_component_interpretation), - postprocessor(nullptr, typeid(*this).name()), - n_output_variables(names.size()) - { - Assert (names.size() == data_component_interpretation.size(), - ExcDimensionMismatch(data_component_interpretation.size(), - names.size())); - - // check that the names use only allowed characters - for (unsigned int i=0; i()") == std::string::npos, - Exceptions::DataOut::ExcInvalidCharacter (names[i], - names[i].find_first_not_of("abcdefghijklmnopqrstuvwxyz" - "ABCDEFGHIJKLMNOPQRSTUVWXYZ" - "0123456789_<>()"))); - } - - - - template - DataEntryBase::DataEntryBase - (const DoFHandlerType *dofs, - const DataPostprocessor *data_postprocessor) - : - dof_handler (dofs, typeid(dealii::DataOut_DoFData).name()), - names(data_postprocessor->get_names()), - data_component_interpretation (data_postprocessor->get_data_component_interpretation()), - postprocessor(data_postprocessor, typeid(*this).name()), - n_output_variables(names.size()) - { - Assert (data_postprocessor->get_names().size() - == - data_postprocessor->get_data_component_interpretation().size(), - ExcDimensionMismatch (data_postprocessor->get_names().size(), - data_postprocessor->get_data_component_interpretation().size())); - - // check that the names use only allowed characters - for (unsigned int i=0; i()") == std::string::npos, - Exceptions::DataOut::ExcInvalidCharacter (names[i], - names[i].find_first_not_of("abcdefghijklmnopqrstuvwxyz" - "ABCDEFGHIJKLMNOPQRSTUVWXYZ" - "0123456789_<>()"))); - } - - - - template - DataEntryBase::~DataEntryBase () - {} - - - - /** - * Class that stores a pointer to a vector of type equal to the template - * argument, and provides the functions to extract data from it. - * - * @author Wolfgang Bangerth, 2004 - */ - template - class DataEntry : public DataEntryBase - { - public: - /** - * Constructor. Give a list of names for the individual components of - * the vector and their interpretation as scalar or vector data. This - * constructor assumes that no postprocessor is going to be used. - */ - DataEntry - (const DoFHandlerType *dofs, - const VectorType *data, - const std::vector &names, - const std::vector &data_component_interpretation); - - /** - * Constructor when a data postprocessor is going to be used. In that - * case, the names and vector declarations are going to be acquired from - * the postprocessor. - */ - DataEntry (const DoFHandlerType *dofs, - const VectorType *data, - const DataPostprocessor *data_postprocessor); - - /** - * Assuming that the stored vector is a cell vector, extract the given - * element from it. - */ - virtual - double - get_cell_data_value (const unsigned int cell_number) const; - - /** - * Given a FEValuesBase object, extract the values on the present cell - * from the vector we actually store. - */ - virtual - void - get_function_values - (const FEValuesBase &fe_patch_values, - std::vector &patch_values) const; - - /** - * Given a FEValuesBase object, extract the values 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_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_hessians - (const FEValuesBase &fe_patch_values, - std::vector > &patch_hessians) 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_hessians - (const FEValuesBase &fe_patch_values, - std::vector > > &patch_hessians_system) const; - - /** - * Clear all references to the vectors. - */ - virtual void clear (); - - /** - * Determine an estimate for the memory consumption (in bytes) of this - * object. - */ - virtual std::size_t memory_consumption () const; - - private: - /** - * Pointer to the data vector. Note that ownership of the vector pointed - * to remains with the caller of this class. - */ - const VectorType *vector; - }; - - - - template - DataEntry:: - DataEntry (const DoFHandlerType *dofs, - const VectorType *data, - const std::vector &names, - const std::vector &data_component_interpretation) - : - DataEntryBase (dofs, names, data_component_interpretation), - vector (data) - {} - - - - template - DataEntry:: - DataEntry (const DoFHandlerType *dofs, - const VectorType *data, - const DataPostprocessor *data_postprocessor) - : - DataEntryBase (dofs, data_postprocessor), - vector (data) - {} - - - namespace - { - template - double - get_vector_element (const VectorType &vector, - const unsigned int cell_number) - { - return internal::ElementAccess::get(vector,cell_number); - } - - - double - get_vector_element (const IndexSet &is, - const unsigned int cell_number) - { - return (is.is_element(cell_number) ? 1 : 0); - } - } - - - - template - double - DataEntry:: - get_cell_data_value (const unsigned int cell_number) const - { - return get_vector_element(*vector, cell_number); - } - - - - template - void - DataEntry::get_function_values - (const FEValuesBase &fe_patch_values, - std::vector > &patch_values_system) const - { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers - if (typeid(typename VectorType::value_type) == typeid(double)) - { - fe_patch_values.get_function_values (*vector, - // reinterpret output argument type; because of - // the 'if' statement above, this is the - // identity cast whenever the code is - // executed, but the cast is necessary - // to allow compilation even if we don't get here - reinterpret_cast >&> - (patch_values_system)); - } - else - { - std::vector > tmp(patch_values_system.size()); - for (unsigned int i = 0; i < patch_values_system.size(); i++) - tmp[i].reinit(patch_values_system[i]); - - fe_patch_values.get_function_values (*vector, tmp); - - for (unsigned int i = 0; i < patch_values_system.size(); i++) - patch_values_system[i] = tmp[i]; - } - } - - - - template - void - DataEntry::get_function_values - (const FEValuesBase &fe_patch_values, - std::vector &patch_values) const - { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers - if (typeid(typename VectorType::value_type) == typeid(double)) - { - fe_patch_values.get_function_values (*vector, - // reinterpret output argument type; because of - // the 'if' statement above, this is the - // identity cast whenever the code is - // executed, but the cast is necessary - // to allow compilation even if we don't get here - reinterpret_cast&> - (patch_values)); - } - else - { - std::vector tmp (patch_values.size()); - - fe_patch_values.get_function_values (*vector, tmp); - - for (unsigned int i = 0; i < tmp.size(); i++) - patch_values[i] = tmp[i]; - } - } - - - - template - void - DataEntry::get_function_gradients - (const FEValuesBase &fe_patch_values, - std::vector > > &patch_gradients_system) const - { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers - if (typeid(typename VectorType::value_type) == typeid(double)) - { - fe_patch_values.get_function_gradients (*vector, - // reinterpret output argument type; because of - // the 'if' statement above, this is the - // identity cast whenever the code is - // executed, but the cast is necessary - // to allow compilation even if we don't get here - reinterpret_cast > >&> - (patch_gradients_system)); - } - else - { - std::vector > > - tmp(patch_gradients_system.size()); - for (unsigned int i = 0; i < tmp.size(); i++) - tmp[i].resize(patch_gradients_system[i].size()); - - fe_patch_values.get_function_gradients (*vector, tmp); - - for (unsigned int i = 0; i < tmp.size(); i++) - for (unsigned int j = 0; j < tmp[i].size(); j++) - patch_gradients_system[i][j] = tmp[i][j]; - } - } - - - - template - void - DataEntry::get_function_gradients - (const FEValuesBase &fe_patch_values, - std::vector > &patch_gradients) const - { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers - if (typeid(typename VectorType::value_type) == typeid(double)) - { - fe_patch_values.get_function_gradients (*vector, - // reinterpret output argument type; because of - // the 'if' statement above, this is the - // identity cast whenever the code is - // executed, but the cast is necessary - // to allow compilation even if we don't get here - reinterpret_cast >&> - (patch_gradients)); - } - else - { - std::vector > tmp; - tmp.resize(patch_gradients.size()); - - fe_patch_values.get_function_gradients (*vector, tmp); - - for (unsigned int i = 0; i < tmp.size(); i++) - patch_gradients[i] = tmp[i]; - } - } - - - - template - void - DataEntry::get_function_hessians - (const FEValuesBase &fe_patch_values, - std::vector > > &patch_hessians_system) const - { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers - if (typeid(typename VectorType::value_type) == typeid(double)) - { - fe_patch_values.get_function_hessians (*vector, - // reinterpret output argument type; because of - // the 'if' statement above, this is the - // identity cast whenever the code is - // executed, but the cast is necessary - // to allow compilation even if we don't get here - reinterpret_cast > >&> - (patch_hessians_system)); - } - else - { - std::vector > > - tmp(patch_hessians_system.size()); - for (unsigned int i = 0; i < tmp.size(); i++) - tmp[i].resize(patch_hessians_system[i].size()); - - fe_patch_values.get_function_hessians (*vector, tmp); - - for (unsigned int i = 0; i < tmp.size(); i++) - for (unsigned int j = 0; j < tmp[i].size(); j++) - patch_hessians_system[i][j] = tmp[i][j]; - } - } - - - - template - void - DataEntry::get_function_hessians - (const FEValuesBase &fe_patch_values, - std::vector > &patch_hessians) const - { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers - if (typeid(typename VectorType::value_type) == typeid(double)) - { - fe_patch_values.get_function_hessians (*vector, - // reinterpret output argument type; because of - // the 'if' statement above, this is the - // identity cast whenever the code is - // executed, but the cast is necessary - // to allow compilation even if we don't get here - reinterpret_cast >&> - (patch_hessians)); - } - else - { - std::vector > - tmp(patch_hessians.size()); - - fe_patch_values.get_function_hessians (*vector, tmp); - - for (unsigned int i = 0; i < tmp.size(); i++) - patch_hessians[i] = tmp[i]; - } - } - - - - template - std::size_t - DataEntry::memory_consumption () const - { - return (sizeof (vector) + - MemoryConsumption::memory_consumption (this->names)); - } - - - - template - void - DataEntry::clear () - { - vector = nullptr; - this->dof_handler = nullptr; - } - } -} - - - -template -DataOut_DoFData::DataOut_DoFData () - : - triangulation(nullptr,typeid(*this).name()), - dofs(nullptr,typeid(*this).name()) -{} - - - -template -DataOut_DoFData::~DataOut_DoFData () -{ - clear (); -} - - - -template -void -DataOut_DoFData:: -attach_dof_handler (const DoFHandlerType &d) -{ - Assert (dof_data.size() == 0, - Exceptions::DataOut::ExcOldDataStillPresent()); - Assert (cell_data.size() == 0, - Exceptions::DataOut::ExcOldDataStillPresent()); - - triangulation = SmartPointer > - (&d.get_triangulation(), typeid(*this).name()); - dofs = SmartPointer(&d, typeid(*this).name()); -} - - - -template -void -DataOut_DoFData:: -attach_triangulation (const Triangulation &tria) -{ - Assert (dof_data.size() == 0, - Exceptions::DataOut::ExcOldDataStillPresent()); - Assert (cell_data.size() == 0, - Exceptions::DataOut::ExcOldDataStillPresent()); - - triangulation = SmartPointer > - (&tria, typeid(*this).name()); -} - - - -template -template -void -DataOut_DoFData:: -add_data_vector (const DoFHandlerType &dof_handler, - const VectorType &vec, - const DataPostprocessor &data_postprocessor) -{ - // this is a specialized version of the other function where we have a - // postprocessor. if we do, we know that we have type_dof_data, which makes - // things a bit simpler, we also don't need to deal with some of the other - // stuff and use a different constructor of DataEntry - if (triangulation != nullptr) - { - Assert (&dof_handler.get_triangulation() == triangulation, - ExcMessage("The triangulation attached to the DoFHandler does not " - "match with the one set previously")); - } - else - { - triangulation = SmartPointer> - (&dof_handler.get_triangulation(), typeid(*this).name()); - } - - Assert (vec.size() == dof_handler.n_dofs(), - Exceptions::DataOut::ExcInvalidVectorSize (vec.size(), - dof_handler.n_dofs(), - dof_handler.get_triangulation().n_active_cells())); - - - auto new_entry = std_cxx14::make_unique> - (&dof_handler, &vec, &data_postprocessor); - dof_data.emplace_back (std::move(new_entry)); -} - - - -template -template -void -DataOut_DoFData:: -add_data_vector_internal -(const DoFHandlerType *dof_handler, - const VectorType &data_vector, - const std::vector &names, - const DataVectorType type, - const std::vector &data_component_interpretation_, - const bool deduce_output_names) -{ - // Check available mesh information: - if (triangulation == nullptr) - { - Assert(dof_handler != nullptr, ExcInternalError()); - triangulation = SmartPointer> - (&dof_handler->get_triangulation(), typeid(*this).name()); - } - - if (dof_handler != nullptr) - { - Assert (&dof_handler->get_triangulation() == triangulation, - ExcMessage("The triangulation attached to the DoFHandler does not " - "match with the one set previously")); - } - - // Figure out the data type: - DataVectorType actual_type = type; - if (type == type_automatic) - { - Assert((dof_handler == nullptr) || (triangulation->n_active_cells() != dof_handler->n_dofs()), - ExcMessage("Unable to determine the type of vector automatically because the number of DoFs " - "is equal to the number of cells. Please specify DataVectorType.")); - - if (data_vector.size() == triangulation->n_active_cells()) - actual_type = type_cell_data; - else - actual_type = type_dof_data; - } - Assert(actual_type == type_cell_data || actual_type == type_dof_data, - ExcInternalError()); - - // If necessary, append '_1', '_2', etc. to component names: - std::vector deduced_names; - if (deduce_output_names && actual_type == type_dof_data) - { - Assert(names.size() == 1, ExcInternalError()); - Assert(dof_handler != nullptr, ExcInternalError()); - Assert(dof_handler->n_dofs() > 0, - ExcMessage("The DoF handler attached to the current output vector " - "does not have any degrees of freedom, so it is not " - "possible to output DoF data in this context.")); - const std::string name = names[0]; - const unsigned int n_components = dof_handler->get_fe().n_components(); - deduced_names.resize (n_components); - if (n_components > 1) - { - for (unsigned int i=0; in_active_cells(), - ExcDimensionMismatch (data_vector.size(), - triangulation->n_active_cells())); - Assert (deduced_names.size() == 1, - Exceptions::DataOut::ExcInvalidNumberOfNames (deduced_names.size(), 1)); - break; - case type_dof_data: - Assert (dof_handler != nullptr, - Exceptions::DataOut::ExcNoDoFHandlerSelected ()); - Assert (data_vector.size() == dof_handler->n_dofs(), - Exceptions::DataOut::ExcInvalidVectorSize (data_vector.size(), - dof_handler->n_dofs(), - triangulation->n_active_cells())); - Assert (deduced_names.size() == dof_handler->get_fe().n_components(), - Exceptions::DataOut::ExcInvalidNumberOfNames (deduced_names.size(), - dof_handler->get_fe().n_components())); - break; - default: - Assert (false, ExcInternalError()); - } - - const auto &data_component_interpretation - = (data_component_interpretation_.size() != 0 - ? - data_component_interpretation_ - : - std::vector - (deduced_names.size(), DataComponentInterpretation::component_is_scalar)); - - // finally, add the data vector: - auto new_entry = std_cxx14::make_unique> - (dof_handler, &data_vector, deduced_names, data_component_interpretation); - - if (actual_type == type_dof_data) - dof_data.emplace_back (std::move(new_entry)); - else - cell_data.emplace_back (std::move(new_entry)); -} - - - -template -void DataOut_DoFData::clear_data_vectors () -{ - dof_data.erase (dof_data.begin(), dof_data.end()); - cell_data.erase (cell_data.begin(), cell_data.end()); - - // delete patches - std::vector dummy; - patches.swap (dummy); -} - - - -template -void -DataOut_DoFData:: -clear_input_data_references () -{ - for (unsigned int i=0; iclear (); - - for (unsigned int i=0; iclear (); - - if (dofs != nullptr) - dofs = nullptr; -} - - - -template -void -DataOut_DoFData::clear () -{ - dof_data.erase (dof_data.begin(), dof_data.end()); - cell_data.erase (cell_data.begin(), cell_data.end()); - - if (dofs != nullptr) - dofs = nullptr; - - // delete patches - std::vector dummy; - patches.swap (dummy); -} - - - -template -std::vector -DataOut_DoFData:: -get_dataset_names () const -{ - std::vector names; - // collect the names of dof - // and cell data - typedef - typename std::vector > >::const_iterator - data_iterator; - - for (data_iterator d=dof_data.begin(); - d!=dof_data.end(); ++d) - for (unsigned int i=0; i<(*d)->names.size(); ++i) - names.push_back ((*d)->names[i]); - for (data_iterator d=cell_data.begin(); d!=cell_data.end(); ++d) - { - Assert ((*d)->names.size() == 1, ExcInternalError()); - names.push_back ((*d)->names[0]); - } - - return names; -} - - - -template -std::vector > -DataOut_DoFData::get_vector_data_ranges () const -{ - std::vector > - ranges; - - // collect the ranges of dof - // and cell data - typedef - typename std::vector > >::const_iterator - data_iterator; - - unsigned int output_component = 0; - for (data_iterator d=dof_data.begin(); - d!=dof_data.end(); ++d) - for (unsigned int i=0; i<(*d)->n_output_variables; - ++i, ++output_component) - // see what kind of data we have - // here. note that for the purpose of - // the current function all we care - // about is vector data - if ((*d)->data_component_interpretation[i] == - DataComponentInterpretation::component_is_part_of_vector) - { - // ensure that there is a - // continuous number of next - // space_dim components that all - // deal with vectors - Assert (i+patch_space_dim <= - (*d)->n_output_variables, - Exceptions::DataOut::ExcInvalidVectorDeclaration (i, - (*d)->names[i])); - for (unsigned int dd=1; dddata_component_interpretation[i+dd] - == - DataComponentInterpretation::component_is_part_of_vector, - Exceptions::DataOut::ExcInvalidVectorDeclaration (i, - (*d)->names[i])); - - // all seems alright, so figure out - // whether there is a common name - // to these components. if not, - // leave the name empty and let the - // output format writer decide what - // to do here - std::string name = (*d)->names[i]; - for (unsigned int dd=1; ddnames[i+dd]) - { - name = ""; - break; - } - - // finally add a corresponding - // range - std::tuple - range (output_component, - output_component+patch_space_dim-1, - name); - - ranges.push_back (range); - - // increase the 'component' counter - // by the appropriate amount, same - // for 'i', since we have already - // dealt with all these components - output_component += patch_space_dim-1; - i += patch_space_dim-1; - } - - // note that we do not have to traverse the - // list of cell data here because cell data - // is one value per (logical) cell and - // therefore cannot be a vector - - // as a final check, the 'component' - // counter should be at the total number of - // components added up now -#ifdef DEBUG - unsigned int n_output_components = 0; - for (data_iterator d=dof_data.begin(); - d!=dof_data.end(); ++d) - n_output_components += (*d)->n_output_variables; - Assert (output_component == n_output_components, - ExcInternalError()); -#endif - - return ranges; -} - - - -template -const std::vector< dealii::DataOutBase::Patch > & -DataOut_DoFData::get_patches () const -{ - return patches; -} - - - -template -std::vector > > - DataOut_DoFData::get_finite_elements() const -{ - const unsigned int dhdim = DoFHandlerType::dimension; - const unsigned int dhspacedim = DoFHandlerType::space_dimension; - std::vector > > - finite_elements(this->dof_data.size()); - for (unsigned int i=0; idof_data.size(); ++i) - { - Assert (dof_data[i]->dof_handler != nullptr, - Exceptions::DataOut::ExcNoDoFHandlerSelected ()); - - // avoid creating too many finite elements and doing a lot of work on - // initializing FEValues downstream: if two DoFHandlers are the same - // (checked by pointer comparison), we can re-use the shared_ptr object - // for the second one. We cannot check for finite element equalities - // because we need different FEValues objects for different dof - // handlers. - bool duplicate = false; - for (unsigned int j=0; jdof_handler == dof_data[j]->dof_handler) - { - finite_elements[i] = finite_elements[j]; - duplicate = true; - } - if (duplicate == false) - finite_elements[i].reset(new dealii::hp::FECollection - (this->dof_data[i]->dof_handler->get_fe())); - } - if (this->dof_data.empty()) - { - finite_elements.resize(1); - finite_elements[0].reset(new dealii::hp::FECollection - (FE_DGQ(0))); - } - return finite_elements; -} - - - -template -std::size_t -DataOut_DoFData::memory_consumption () const -{ - return (DataOutInterface::memory_consumption () + - MemoryConsumption::memory_consumption (dofs) + - MemoryConsumption::memory_consumption (patches)); -} - - -// explicit instantiations #include "data_out_dof_data.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/source/numerics/data_out_dof_data.inst.in b/source/numerics/data_out_dof_data.inst.in index 32bf5f57fd..8ceb65aacb 100644 --- a/source/numerics/data_out_dof_data.inst.in +++ b/source/numerics/data_out_dof_data.inst.in @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2010 - 2016 by the deal.II authors +// Copyright (C) 2010 - 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -16,8 +16,6 @@ for (VEC : REAL_SERIAL_VECTORS; DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) { -// codim=0 - template void DataOut_DoFData,deal_II_dimension,deal_II_dimension>:: add_data_vector_internal (const DH *, @@ -72,56 +70,12 @@ for (VEC : REAL_SERIAL_VECTORS; DH : DOFHANDLER_TEMPLATES; deal_II_dimension : D const VEC &, const DataPostprocessor::space_dimension> &); #endif - - -// codim 1 - -#if deal_II_dimension < 3 - template void - DataOut_DoFData,deal_II_dimension,deal_II_dimension+1>:: - add_data_vector_internal (const DH *, - const VEC &, - const std::vector &, - const DataVectorType, - const std::vector &, - const bool); - - template void - DataOut_DoFData,deal_II_dimension,deal_II_dimension+1>:: - add_data_vector (const DH &, - const VEC &, - const DataPostprocessor::space_dimension> &); -#endif - - - -// codim 2 - -#if deal_II_dimension == 3 - template void - DataOut_DoFData,1,3>:: - add_data_vector_internal (const DH<1, 3> *, - const VEC &, - const std::vector &, - const DataVectorType, - const std::vector &, - const bool); - - template void - DataOut_DoFData,1,3>:: - add_data_vector (const DH<1,3> &, - const VEC &, - const DataPostprocessor::space_dimension> &); -#endif - } for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) { -// codim=0 - template void DataOut_DoFData,deal_II_dimension,deal_II_dimension>:: add_data_vector_internal (const DH *, @@ -133,9 +87,9 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) template void DataOut_DoFData,deal_II_dimension,deal_II_dimension>:: - add_data_vector (const DH &, - const IndexSet &, - const DataPostprocessor::space_dimension> &); + add_data_vector (const DH &, + const IndexSet &, + const DataPostprocessor::space_dimension> &); @@ -174,26 +128,6 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) const IndexSet &, const DataPostprocessor::space_dimension> &); #endif - -// codim 1 - -#if deal_II_dimension < 3 - template void - DataOut_DoFData,deal_II_dimension,deal_II_dimension+1>:: - add_data_vector_internal (const DH *, - const IndexSet &, - const std::vector &, - const DataVectorType, - const std::vector &, - const bool); - - template void - DataOut_DoFData,deal_II_dimension,deal_II_dimension+1>:: - add_data_vector (const DH &, - const IndexSet &, - const DataPostprocessor::space_dimension> &); -#endif - } @@ -205,7 +139,6 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) #if deal_II_dimension < 3 template class DataOut_DoFData,deal_II_dimension+1>; template class DataOut_DoFData,deal_II_dimension,deal_II_dimension+1>; - template class DataOut_DoFData,deal_II_dimension,deal_II_dimension+1>; #endif #if deal_II_dimension >= 2 @@ -214,37 +147,32 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) #if deal_II_dimension == 3 template class DataOut_DoFData,1,3>; - template class DataOut_DoFData,1,3>; #endif } -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) +for (deal_II_dimension : DIMENSIONS) { namespace internal \{ namespace DataOut \{ -#if deal_II_dimension <= deal_II_space_dimension - template struct ParallelDataBase; -#endif + template struct ParallelDataBase; \} \} } -for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) +for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) { namespace internal \{ namespace DataOut \{ -#if deal_II_dimension <= deal_II_space_dimension template void - ParallelDataBase:: - reinit_all_fe_values > - (std::vector > > > &dof_data, - const dealii::Triangulation::cell_iterator &cell, + ParallelDataBase:: + reinit_all_fe_values > + (std::vector > > > &dof_data, + const dealii::Triangulation::cell_iterator &cell, const unsigned int face); -#endif \} \} } diff --git a/source/numerics/data_out_dof_data_codim.cc b/source/numerics/data_out_dof_data_codim.cc new file mode 100644 index 0000000000..bbb244d16c --- /dev/null +++ b/source/numerics/data_out_dof_data_codim.cc @@ -0,0 +1,22 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include + +DEAL_II_NAMESPACE_OPEN + +#include "data_out_dof_data_codim.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/numerics/data_out_dof_data_codim.inst.in b/source/numerics/data_out_dof_data_codim.inst.in new file mode 100644 index 0000000000..f7e709b7af --- /dev/null +++ b/source/numerics/data_out_dof_data_codim.inst.in @@ -0,0 +1,106 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +for (VEC : REAL_SERIAL_VECTORS; DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) +{ +#if deal_II_dimension < deal_II_space_dimension + template void + DataOut_DoFData,deal_II_dimension,deal_II_space_dimension>:: + add_data_vector_internal (const DH *, + const VEC &, + const std::vector &, + const DataVectorType, + const std::vector &, + const bool); + + template void + DataOut_DoFData,deal_II_dimension,deal_II_space_dimension>:: + add_data_vector (const VEC &, + const DataPostprocessor::space_dimension> &); + + template void + DataOut_DoFData,deal_II_dimension,deal_II_space_dimension>:: + add_data_vector (const DH &, + const VEC &, + const DataPostprocessor::space_dimension> &); +#endif +} + + + +for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) +{ +#if deal_II_dimension < deal_II_space_dimension + template void + DataOut_DoFData,deal_II_dimension,deal_II_space_dimension>:: + add_data_vector_internal + (const DH *, + const IndexSet &, + const std::vector &, + const DataVectorType, + const std::vector &, + const bool); + + template void + DataOut_DoFData,deal_II_dimension,deal_II_space_dimension>:: + add_data_vector (const IndexSet &, + const DataPostprocessor::space_dimension> &); +#endif +} + + + +for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) +{ +#if deal_II_dimension < 3 + template class DataOut_DoFData,deal_II_dimension,deal_II_dimension+1>; +#endif + +#if deal_II_dimension == 3 + template class DataOut_DoFData,1,3>; +#endif + +} + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) +{ + namespace internal \{ + namespace DataOut \{ +#if deal_II_dimension < deal_II_space_dimension + template struct ParallelDataBase; +#endif + \} + \} +} + + +for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) +{ + namespace internal \{ + namespace DataOut \{ +#if deal_II_dimension < deal_II_space_dimension + template + void + ParallelDataBase:: + reinit_all_fe_values > + (std::vector > > > &dof_data, + const dealii::Triangulation::cell_iterator &cell, + const unsigned int face); +#endif + \} + \} +} -- 2.39.5