From: Peter Munch Date: Thu, 4 Mar 2021 16:37:28 +0000 (+0100) Subject: Use QCollection within MapptingInfo::initialize X-Git-Tag: v9.3.0-rc1~353^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=83e5263735583d42fe4227e0dd23630f8e0dd0b1;p=dealii.git Use QCollection within MapptingInfo::initialize --- diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index 8c0d0c77d4..4ba8c6b73c 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -27,6 +27,8 @@ #include #include +#include + #include #include @@ -194,6 +196,14 @@ namespace internal */ std::vector descriptor; + /** + * Collection of quadrature formulae applied on the given face. + * + * @note Only filled for faces, since faces might be quadrilateral or + * triangle shaped. + */ + std::vector> q_collection; + /** * Stores the index offset into the arrays @p jxw_values, @p jacobians, * @p normal_vectors and the second derivatives. Note that affine cells @@ -454,6 +464,12 @@ namespace internal */ SmartPointer> mapping; + /** + * Reference-cell type related to each quadrature and active quadrature + * index. + */ + std::vector> reference_cell_types; + /** * Internal function to compute the geometry for the case the mapping is * a MappingQ and a single quadrature formula per slot (non-hp-case) is diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 445e28dea7..d9d2f84522 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -362,37 +362,21 @@ namespace internal this->update_flags_inner_faces = this->update_flags_boundary_faces; this->update_flags_faces_by_cells = update_flags_faces_by_cells; + reference_cell_types.resize(quad.size()); + for (unsigned int my_q = 0; my_q < quad.size(); ++my_q) { const unsigned int n_hp_quads = quad[my_q].size(); AssertIndexRange(0, n_hp_quads); - cell_data[my_q].descriptor.resize(n_hp_quads); - for (unsigned int q = 0; q < n_hp_quads; ++q) - { - bool flag = quad[my_q][q].is_tensor_product(); - if (flag) - for (unsigned int i = 1; i < dim; ++i) - { - flag &= quad[my_q][q].get_tensor_basis()[0] == - quad[my_q][q].get_tensor_basis()[i]; - } + const unsigned int scale = std::max(1, dim - 1); - if (flag == false) - { -#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT - cell_data[my_q].descriptor[q].initialize(quad[my_q][q], - update_default); -#else - Assert(false, ExcNotImplemented()); -#endif - } - else - cell_data[my_q].descriptor[q].initialize( - quad[my_q][q].get_tensor_basis()[0], update_default); - } + cell_data[my_q].descriptor.resize(n_hp_quads); + face_data[my_q].descriptor.resize(n_hp_quads * scale); + face_data[my_q].q_collection.resize(n_hp_quads); + face_data_by_cells[my_q].descriptor.resize(n_hp_quads * scale); + reference_cell_types[my_q].resize(n_hp_quads); - face_data[my_q].descriptor.resize(n_hp_quads); for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq) { bool flag = quad[my_q][hpq].is_tensor_product(); @@ -405,58 +389,57 @@ namespace internal if (flag == false) { #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT - try + cell_data[my_q].descriptor[hpq].initialize(quad[my_q][hpq], + update_default); + const auto quad_face = + get_unique_face_quadratures(quad[my_q][hpq]); + + if (quad_face.first.size() > 0) // line, quad { - const auto quad_face = - get_face_quadrature(quad[my_q][hpq]); - face_data[my_q].descriptor[hpq].initialize( - quad_face, update_default); + face_data[my_q].descriptor[hpq * scale].initialize( + quad_face.first, update_default); + face_data_by_cells[my_q] + .descriptor[hpq * scale] + .initialize(quad_face.first, update_default); } - catch (...) + + if (quad_face.second.size() > 0) // triangle { - // TODO: nothing to do for now for wedges and pyramids. + AssertDimension(dim, 3); + face_data[my_q].descriptor[hpq * scale + 1].initialize( + quad_face.second, update_default); + face_data_by_cells[my_q] + .descriptor[hpq * scale + 1] + .initialize(quad_face.second, update_default); } -#else - Assert(false, ExcNotImplemented()); -#endif - } - else - face_data[my_q].descriptor[hpq].initialize( - quad[my_q][hpq].get_tensor_basis()[0], - update_flags_boundary_faces); - } - face_data_by_cells[my_q].descriptor.resize(n_hp_quads); - for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq) - { - bool flag = quad[my_q][hpq].is_tensor_product(); + const auto face_quadrature_collection = + get_face_quadrature_collection(quad[my_q][hpq]); - if (flag) - for (unsigned int i = 1; i < dim; ++i) - flag &= quad[my_q][hpq].get_tensor_basis()[0] == - quad[my_q][hpq].get_tensor_basis()[i]; + face_data[my_q].q_collection[hpq] = + face_quadrature_collection.second; - if (flag == false) - { -#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT - try - { - const auto quad_face = - get_face_quadrature(quad[my_q][hpq]); - face_data_by_cells[my_q].descriptor[hpq].initialize( - quad_face, update_default); - } - catch (...) - { - // TODO: nothing to do for now for wedges and pyramids. - } + reference_cell_types[my_q][hpq] = + face_quadrature_collection.first; #else Assert(false, ExcNotImplemented()); #endif } else - face_data_by_cells[my_q].descriptor[hpq].initialize( - quad[my_q][hpq].get_tensor_basis()[0], update_default); + { + cell_data[my_q].descriptor[hpq].initialize( + quad[my_q][hpq].get_tensor_basis()[0], update_default); + + const auto quad_face = quad[my_q][hpq].get_tensor_basis()[0]; + face_data[my_q].descriptor[hpq * scale].initialize( + quad_face, update_flags_boundary_faces); + face_data[my_q].q_collection[hpq] = + dealii::hp::QCollection(quad_face); + face_data_by_cells[my_q].descriptor[hpq * scale].initialize( + quad_face, update_default); + reference_cell_types[my_q][hpq] = + dealii::ReferenceCells::get_hypercube(); + } } } @@ -1804,7 +1787,20 @@ namespace internal MappingInfoStorage>, CompressedFaceData> &data) { - FE_Nothing dummy_fe; + std::vector>>> dummy_fe( + mapping_info.reference_cell_types.size()); + for (unsigned int my_q = 0; + my_q < mapping_info.reference_cell_types.size(); + ++my_q) + { + const unsigned int n_hp_quads = + mapping_info.reference_cell_types[my_q].size(); + dummy_fe[my_q].resize(n_hp_quads); + + for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq) + dummy_fe[my_q][hpq] = std::make_shared>( + mapping_info.reference_cell_types[my_q][hpq]); + } const unsigned int max_active_fe_index = active_fe_index.size() > 0 ? @@ -1860,11 +1856,9 @@ namespace internal const unsigned int hp_mapping_index = mapping_in.size() == 1 ? 0 : fe_index; - const auto & mapping = mapping_in[hp_mapping_index]; - const Quadrature &quadrature = - mapping_info.face_data[my_q] - .descriptor[hp_quad_index] - .quadrature; + const auto &mapping = mapping_in[hp_mapping_index]; + const auto &quadrature = + mapping_info.face_data[my_q].q_collection[hp_quad_index]; const bool is_boundary_face = faces[face].cells_exterior[0] == numbers::invalid_unsigned_int; @@ -1874,14 +1868,14 @@ namespace internal fe_boundary_face_values_container[my_q][fe_index] = std::make_shared>( mapping, - dummy_fe, + *dummy_fe[my_q][hp_quad_index], quadrature, mapping_info.update_flags_boundary_faces); else if (fe_face_values_container[my_q][fe_index] == nullptr) fe_face_values_container[my_q][fe_index] = std::make_shared>( mapping, - dummy_fe, + *dummy_fe[my_q][hp_quad_index], quadrature, mapping_info.update_flags_inner_faces); @@ -1889,9 +1883,9 @@ namespace internal is_boundary_face ? *fe_boundary_face_values_container[my_q][fe_index] : *fe_face_values_container[my_q][fe_index]; - const unsigned int n_q_points = - fe_face_values.n_quadrature_points; - face_data.resize(n_q_points); + + unsigned int n_q_points = 0; // will be override once FEFaceValues + // is set up bool normal_is_similar = true; bool JxW_is_similar = true; @@ -1911,13 +1905,22 @@ namespace internal fe_face_values.reinit(cell_it, faces[face].interior_face_no); + if (v == 0) + { + n_q_points = fe_face_values.n_quadrature_points; + face_data.resize(n_q_points); + } + for (unsigned int q = 0; q < n_q_points; ++q) { if (std::abs( - fe_face_values.JxW(q) * quadrature.weight(0) - - fe_face_values.JxW(0) * quadrature.weight(q)) > + fe_face_values.JxW(q) * + fe_face_values.get_quadrature().weight(0) - + fe_face_values.JxW(0) * + fe_face_values.get_quadrature().weight(q)) > 2048. * std::numeric_limits::epsilon() * - fe_face_values.JxW(0) * quadrature.weight(q)) + fe_face_values.JxW(0) * + fe_face_values.get_quadrature().weight(q)) JxW_is_similar = false; face_data.JxW_values[q][v] = fe_face_values.JxW(q); @@ -2043,17 +2046,16 @@ namespace internal mapping_in.size() == 1 ? 0 : fe_index; const auto &mapping = mapping_in[hp_mapping_index]; - const Quadrature &quadrature = + const auto &quadrature = mapping_info.face_data[my_q] - .descriptor[hp_quad_index] - .quadrature; + .q_collection[hp_quad_index]; if (fe_face_values_container[my_q][fe_index] == nullptr) fe_face_values_container[my_q][fe_index] = std::make_shared>( mapping, - dummy_fe, + *dummy_fe[my_q][hp_quad_index], quadrature, mapping_info.update_flags_boundary_faces); @@ -2069,7 +2071,7 @@ namespace internal fe_subface_values_container[my_q][0] = std::make_shared>( mapping, - dummy_fe, + *dummy_fe[my_q][hp_quad_index], quadrature, mapping_info.update_flags_inner_faces); fe_subface_values_container[my_q][0]->reinit( @@ -2178,7 +2180,8 @@ namespace internal for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) new_entry.first[2 * dim * dim + dim][v] = - face_data.JxW_values[0][v] / quadrature.weight(0) / + face_data.JxW_values[0][v] / + fe_face_values.get_quadrature().weight(0) / Utilities::fixed_power(face_data.jac_size); new_entry.second = data.second.data.size(); @@ -3165,14 +3168,14 @@ namespace internal std::make_shared>( mapping, dummy_fe, - face_data_by_cells[my_q].descriptor[fe_index].quadrature, + face_data[my_q].q_collection[fe_index], update_flags); if (fe_face_values_neigh[my_q][fe_index].get() == nullptr) fe_face_values_neigh[my_q][fe_index] = std::make_shared>( mapping, dummy_fe, - face_data_by_cells[my_q].descriptor[fe_index].quadrature, + face_data[my_q].q_collection[fe_index], update_flags); dealii::FEFaceValues &fe_val = *fe_face_values[my_q][fe_index]; @@ -3215,9 +3218,7 @@ namespace internal { if (update_flags & update_JxW_values) face_data_by_cells[my_q].JxW_values[offset][v] = - fe_val.JxW(0) / face_data_by_cells[my_q] - .descriptor[fe_index] - .quadrature.weight(0); + fe_val.JxW(0) / fe_val.get_quadrature().weight(0); if (update_flags & update_jacobians) { DerivativeForm<1, dim, dim> inv_jac = diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h index 164b759e27..2c47b96fab 100644 --- a/include/deal.II/matrix_free/util.h +++ b/include/deal.II/matrix_free/util.h @@ -34,20 +34,6 @@ namespace internal { namespace MatrixFreeFunctions { - template - inline Quadrature - get_face_quadrature(const Quadrature &quad) - { - if (dim == 2 || dim == 3) - for (unsigned int i = 1; i <= 3; ++i) - if (quad == QGaussSimplex(i)) - return QGaussSimplex(i); - - AssertThrow(false, ExcNotImplemented()); - - return Quadrature(); - } - template inline std::pair> get_face_quadrature_collection(const Quadrature &quad, @@ -112,6 +98,36 @@ namespace internal return {ReferenceCells::Invalid, dealii::hp::QCollection()}; } + + + template + inline std::pair, Quadrature> + get_unique_face_quadratures(const Quadrature &quad) + { + if (dim == 2 || dim == 3) + for (unsigned int i = 1; i <= 3; ++i) + if (quad == QGaussSimplex(i)) + { + if (dim == 2) + return {QGaussSimplex(i), Quadrature()}; + else + return {Quadrature(), QGaussSimplex(i)}; + } + + if (dim == 3) + for (unsigned int i = 1; i <= 3; ++i) + if (quad == QGaussWedge(i)) + return {QGauss(i), QGaussSimplex(i)}; + + if (dim == 3) + for (unsigned int i = 1; i <= 2; ++i) + if (quad == QGaussPyramid(i)) + return {QGauss(i), QGaussSimplex(i)}; + + AssertThrow(false, ExcNotImplemented()); + + return {QGauss(1), QGauss(1)}; + } } // end of namespace MatrixFreeFunctions } // end of namespace internal diff --git a/tests/simplex/matrix_free_03.cc b/tests/simplex/matrix_free_03.cc index 0cf57e3f63..131750788c 100644 --- a/tests/simplex/matrix_free_03.cc +++ b/tests/simplex/matrix_free_03.cc @@ -169,8 +169,10 @@ public: } }, [&](const auto &data, auto &dst, const auto &src, const auto face_range) { - FEFaceEvaluation fe_eval(data, true); - FEFaceEvaluation fe_eval_neighbor(data, false); + FEFaceEvaluation fe_eval(data, face_range, true); + FEFaceEvaluation fe_eval_neighbor(data, + face_range, + false); for (unsigned int face = face_range.first; face < face_range.second; face++) @@ -214,7 +216,7 @@ public: } }, [&](const auto &data, auto &dst, const auto &src, const auto face_range) { - FEFaceEvaluation fe_eval(data, true); + FEFaceEvaluation fe_eval(data, face_range, true); for (unsigned int face = face_range.first; face < face_range.second; face++) {