From: Peter Munch Date: Sat, 28 Nov 2020 15:10:11 +0000 (+0100) Subject: Enable DataOut for Wedges X-Git-Tag: v9.3.0-rc1~775^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a635393d18737e57a631cde4dbc83e59321832b4;p=dealii.git Enable DataOut for Wedges --- diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h index f7b9b6fbe2..e9f7ee2b61 100644 --- a/include/deal.II/numerics/data_out_dof_data.templates.h +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -82,78 +82,94 @@ namespace internal {} + + /** + * Set up vectors of FEValues and FEFaceValues needed inside of + * ParallelDataBase and return the maximum number of quadrature points + * needed to allocate enough memory for the scratch data. + */ template - ParallelDataBase::ParallelDataBase( - const unsigned int n_datasets, - const unsigned int n_subdivisions, - const std::vector &n_postprocessor_outputs, - const dealii::hp::MappingCollection &mapping, + unsigned int + setup_parallel_data_base_internal( + const unsigned int n_subdivisions, + const dealii::hp::MappingCollection &mapping_collection, const std::vector< std::shared_ptr>> & 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) + const bool use_face_values, + std::vector>> + &x_fe_values, + std::vector>> + &x_fe_face_values) { unsigned int n_q_points = 0; if (use_face_values == false) { - const bool do_simplex = - std::any_of(finite_elements.begin(), - finite_elements.end(), - [](const auto &fe) { - for (unsigned int i = 0; i < fe->size(); ++i) - if ((*fe)[i].reference_cell_type() != - ReferenceCell::get_hypercube(dim)) - return true; - return false; - }); - - const bool do_hex = - std::any_of(finite_elements.begin(), - finite_elements.end(), - [](const auto &fe) { - for (unsigned int i = 0; i < fe->size(); ++i) - if ((*fe)[i].reference_cell_type() == - ReferenceCell::get_hypercube(dim)) - return true; - return false; - }); - - std::unique_ptr> quadrature1; - std::unique_ptr> quadrature2; - - if (do_simplex) + // determine if specific quadrature rules need to set up + bool needs_hypercube_setup = false; + bool needs_simplex_setup = false; + bool needs_wedge_setup = false; + + for (const auto &fe : finite_elements) + for (unsigned int i = 0; i < fe->size(); ++i) + switch ((*fe)[i].reference_cell_type()) + { + case ReferenceCell::Type::Vertex: + case ReferenceCell::Type::Line: + case ReferenceCell::Type::Quad: + case ReferenceCell::Type::Hex: + needs_hypercube_setup |= true; + break; + case ReferenceCell::Type::Tri: + case ReferenceCell::Type::Tet: + needs_simplex_setup |= true; + break; + case ReferenceCell::Type::Wedge: + needs_wedge_setup |= true; + break; + default: + Assert(false, ExcNotImplemented()); + } + + std::unique_ptr> quadrature_simplex; + std::unique_ptr> quadrature_hypercube; + std::unique_ptr> quadrature_wedge; + + if (needs_simplex_setup) { - Assert(1 <= n_subdivisions && n_subdivisions <= 2, - ExcNotImplemented()); - quadrature1.reset(new Quadrature( - Simplex::FE_P(n_subdivisions == 1 ? 1 : 2) - .get_unit_support_points())); + quadrature_simplex = std::make_unique>( + Simplex::FE_P(n_subdivisions) + .get_unit_support_points()); } - if (do_hex) + if (needs_hypercube_setup) { - quadrature2.reset(new Quadrature( - QIterated(QTrapezoid<1>(), n_subdivisions))); + quadrature_hypercube = + std::make_unique>(QTrapezoid<1>(), + n_subdivisions); } - n_q_points = std::max(do_simplex ? quadrature1->size() : 0, - do_hex ? quadrature2->size() : 1); + if (needs_wedge_setup) + { + quadrature_wedge = std::make_unique>( + Simplex::FE_WedgeP( + 1 /*note: vtk only supports linear wedges*/) + .get_unit_support_points()); + } + + n_q_points = std::max( + {needs_wedge_setup ? quadrature_wedge->size() : 0, + needs_simplex_setup ? quadrature_simplex->size() : 0, + needs_hypercube_setup ? quadrature_hypercube->size() : 0}); - x_fe_values.resize(this->finite_elements.size()); - for (unsigned int i = 0; i < this->finite_elements.size(); ++i) + x_fe_values.resize(finite_elements.size()); + for (unsigned int i = 0; i < finite_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 + // 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; j < i; ++j) - if (this->finite_elements[i].get() == - this->finite_elements[j].get()) + if (finite_elements[i].get() == finite_elements[j].get()) { x_fe_values[i] = x_fe_values[j]; break; @@ -162,21 +178,37 @@ namespace internal { dealii::hp::QCollection quadrature; - for (unsigned int j = 0; j < this->finite_elements[i]->size(); - ++j) - if ((*this->finite_elements[i])[j].reference_cell_type() != - ReferenceCell::get_hypercube(dim)) - quadrature.push_back(*quadrature1); - else - quadrature.push_back(*quadrature2); + for (unsigned int j = 0; j < finite_elements[i]->size(); ++j) + switch ((*finite_elements[i])[j].reference_cell_type()) + { + case ReferenceCell::Type::Vertex: + case ReferenceCell::Type::Line: + case ReferenceCell::Type::Quad: + case ReferenceCell::Type::Hex: + quadrature.push_back(*quadrature_hypercube); + break; + case ReferenceCell::Type::Tri: + case ReferenceCell::Type::Tet: + quadrature.push_back(*quadrature_simplex); + break; + case ReferenceCell::Type::Wedge: + quadrature.push_back(*quadrature_wedge); + break; + default: + Assert(false, ExcNotImplemented()); + } x_fe_values[i] = std::make_shared>( - this->mapping_collection, - *this->finite_elements[i], + mapping_collection, + *finite_elements[i], quadrature, - do_simplex ? (update_flags | update_quadrature_points) : - update_flags); + // certain reference-cell kinds take by default the + // curved code path in DataOut::build_one_patch, for + // which the quadrature points need to be enabled + (needs_simplex_setup || needs_wedge_setup) ? + (update_flags | update_quadrature_points) : + update_flags); } } } @@ -185,14 +217,13 @@ namespace internal dealii::hp::QCollection quadrature( QIterated(QTrapezoid<1>(), n_subdivisions)); n_q_points = quadrature[0].size(); - x_fe_face_values.resize(this->finite_elements.size()); - for (unsigned int i = 0; i < this->finite_elements.size(); ++i) + x_fe_face_values.resize(finite_elements.size()); + for (unsigned int i = 0; i < finite_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 + // 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; j < i; ++j) - if (this->finite_elements[i].get() == - this->finite_elements[j].get()) + if (finite_elements[i].get() == finite_elements[j].get()) { x_fe_face_values[i] = x_fe_face_values[j]; break; @@ -200,13 +231,45 @@ namespace internal if (x_fe_face_values[i].get() == nullptr) x_fe_face_values[i] = std::make_shared>( - this->mapping_collection, - *this->finite_elements[i], + mapping_collection, + *finite_elements[i], quadrature, - this->update_flags); + update_flags); } } + return n_q_points; + } + + + + template + ParallelDataBase::ParallelDataBase( + const unsigned int n_datasets, + const unsigned int n_subdivisions, + const std::vector &n_postprocessor_outputs, + const dealii::hp::MappingCollection &mapping, + const std::vector< + std::shared_ptr>> + & 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) + { + const unsigned int n_q_points = + setup_parallel_data_base_internal(n_subdivisions, + mapping, + finite_elements, + update_flags, + use_face_values, + x_fe_values, + x_fe_face_values); + 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); @@ -238,110 +301,13 @@ namespace internal , 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()); - - const bool do_simplex = - std::any_of(finite_elements.begin(), - finite_elements.end(), - [](const auto &fe) { - for (unsigned int i = 0; i < fe->size(); ++i) - if ((*fe)[i].reference_cell_type() != - ReferenceCell::get_hypercube(dim)) - return true; - return false; - }); - - const bool do_hex = - std::any_of(finite_elements.begin(), - finite_elements.end(), - [](const auto &fe) { - for (unsigned int i = 0; i < fe->size(); ++i) - if ((*fe)[i].reference_cell_type() == - ReferenceCell::get_hypercube(dim)) - return true; - return false; - }); - - std::unique_ptr> quadrature1; - std::unique_ptr> quadrature2; - - if (do_simplex) - { - Assert(1 <= n_subdivisions && n_subdivisions <= 2, - ExcNotImplemented()); - quadrature1.reset(new Quadrature( - Simplex::FE_P(n_subdivisions == 1 ? 1 : 2) - .get_unit_support_points())); - } - - if (do_hex) - { - quadrature2.reset(new Quadrature( - QIterated(QTrapezoid<1>(), n_subdivisions))); - } - - x_fe_values.resize(this->finite_elements.size()); - - for (unsigned int i = 0; i < this->finite_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; j < i; ++j) - if (this->finite_elements[i].get() == - this->finite_elements[j].get()) - { - x_fe_values[i] = x_fe_values[j]; - break; - } - if (x_fe_values[i].get() == nullptr) - { - dealii::hp::QCollection quadrature; - - for (unsigned int j = 0; j < this->finite_elements[i]->size(); - ++j) - if ((*this->finite_elements[i])[j].reference_cell_type() != - ReferenceCell::get_hypercube(dim)) - quadrature.push_back(*quadrature1); - else - quadrature.push_back(*quadrature2); - - x_fe_values[i] = - std::make_shared>( - this->mapping_collection, - *this->finite_elements[i], - quadrature, - do_simplex ? (update_flags | update_quadrature_points) : - update_flags); - } - } - } - else - { - dealii::hp::QCollection quadrature( - QIterated(QTrapezoid<1>(), n_subdivisions)); - x_fe_face_values.resize(this->finite_elements.size()); - for (unsigned int i = 0; i < this->finite_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; j < i; ++j) - if (this->finite_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] = - std::make_shared>( - this->mapping_collection, - *this->finite_elements[i], - quadrature, - this->update_flags); - } - } + setup_parallel_data_base_internal(n_subdivisions, + mapping_collection, + finite_elements, + update_flags, + data.x_fe_values.empty(), + x_fe_values, + x_fe_face_values); } diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index 0d91d13b26..f3a51f184b 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -628,6 +628,12 @@ namespace vtk_cell_id[0] = 24; vtk_cell_id[1] = 1; } + else if (patch.reference_cell_type == ReferenceCell::Type::Wedge && + patch.data.n_cols() == 6) + { + vtk_cell_id[0] = 13; + vtk_cell_id[1] = 1; + } else if (patch.reference_cell_type == ReferenceCell::get_hypercube(dim)) { vtk_cell_id[0] = vtk_cell_type[dim];