From: Martin Kronbichler Date: Mon, 13 Dec 2021 11:51:56 +0000 (+0100) Subject: Get rid of std::tuple by providing simple constructor and InitializationData X-Git-Tag: v9.4.0-rc1~740^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=df757529c6aab8af0d163a852914d7ec398f3cd0;p=dealii.git Get rid of std::tuple by providing simple constructor and InitializationData --- diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index c2ad7ed262..3598ec1b60 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -1130,13 +1130,7 @@ namespace internal return; } - internal::MatrixFreeFunctions:: - MappingInfoStorage - temp_data; - temp_data.descriptor.resize(1); - temp_data.descriptor[0].n_q_points = n_q_points; - FEEvaluationData eval( - std::make_tuple(&data.shape_info, nullptr, &temp_data, 0, 0)); + FEEvaluationData eval(data.shape_info); // prepare arrays if (evaluation_flag != EvaluationFlags::nothing) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 1467f754cb..72327cd070 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2972,7 +2972,7 @@ namespace internal dim, n_components, evaluation_flag, - &eval.get_orientation_map()(face_orientation, 0), + &eval.get_shape_info().face_orientations_quad(face_orientation, 0), false, data.n_q_points_face, temp, @@ -3080,7 +3080,7 @@ namespace internal dim, n_components, integration_flag, - &eval.get_orientation_map()(face_orientation, 0), + &eval.get_shape_info().face_orientations_quad(face_orientation, 0), true, data.n_q_points_face, temp, @@ -3203,13 +3203,12 @@ namespace internal if (shape_data.nodal_at_cell_boundaries && eval.get_all_face_orientations()[v] != 0) - orientation[v] = - &eval.get_shape_info() - .face_orientations[eval.get_all_face_orientations()[v]][0]; + orientation[v] = &eval.get_shape_info().face_orientations_dofs( + eval.get_all_face_orientations()[v], 0); } else if (eval.get_face_orientation() != 0) - orientation[0] = - &eval.get_orientation_map()(eval.get_face_orientation(), 0); + orientation[0] = &eval.get_shape_info().face_orientations_dofs( + eval.get_face_orientation(), 0); // face_to_cell_index_hermite std::array index_array_hermite = @@ -3489,7 +3488,8 @@ namespace internal const bool vectorization_possible = all_faces_are_same && (sm_ptr == nullptr); - std::array vector_ptrs; + std::array vector_ptrs; + std::array reordered_indices; if (vectorization_possible == false) { @@ -3543,6 +3543,14 @@ namespace internal Assert(false, ExcNotImplemented()); } } + else if (n_face_orientations == n_lanes) + { + for (unsigned int v = 0; v < n_lanes; ++v) + reordered_indices[v] = + dof_info.dof_indices_contiguous[dof_access_index] + [eval.get_cell_ids()[v]]; + dof_indices = reordered_indices.data(); + } if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients)) { @@ -3741,8 +3749,8 @@ namespace internal n_components, v, evaluation_flag, - &eval.get_orientation_map() - [eval.get_all_face_orientations()[v]][0], + &eval.get_shape_info().face_orientations_quad( + eval.get_all_face_orientations()[v], 0), false, Utilities::pow(n_q_points_1d, dim - 1), &temp[0][0], @@ -3757,7 +3765,8 @@ namespace internal dim, n_components, evaluation_flag, - &eval.get_orientation_map()[eval.get_face_orientation()][0], + &eval.get_shape_info().face_orientations_quad( + eval.get_face_orientation(), 0), false, Utilities::pow(n_q_points_1d, dim - 1), temp, @@ -3918,8 +3927,8 @@ namespace internal n_components, v, integration_flag, - &eval.get_orientation_map() - [eval.get_all_face_orientations()[v]][0], + &eval.get_shape_info().face_orientations_quad( + eval.get_all_face_orientations()[v], 0), true, Utilities::pow(n_q_points_1d, dim - 1), &temp[0][0], @@ -3932,7 +3941,8 @@ namespace internal dim, n_components, integration_flag, - &eval.get_orientation_map()[eval.get_face_orientation()][0], + &eval.get_shape_info().face_orientations_quad( + eval.get_face_orientation(), 0), true, Utilities::pow(n_q_points_1d, dim - 1), temp, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index f05d1b2ac4..2a860a4a8c 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2806,51 +2806,58 @@ namespace internal int dim, typename Number, typename VectorizedArrayType> - inline std::tuple< - const internal::MatrixFreeFunctions::ShapeInfo *, - const internal::MatrixFreeFunctions::DoFInfo *, - const internal::MatrixFreeFunctions:: - MappingInfoStorage *, - unsigned int, - unsigned int> - get_shape_info_and_indices( - const MatrixFree &data, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int fe_degree, - const unsigned int n_q_points, - const unsigned int active_fe_index_given, - const unsigned int active_quad_index_given) + inline typename FEEvaluationData:: + InitializationData + get_shape_info_and_indices( + const MatrixFree &data, + const unsigned int dof_no, + const unsigned int first_selected_component, + const unsigned int quad_no, + const unsigned int fe_degree, + const unsigned int n_q_points, + const unsigned int active_fe_index_given, + const unsigned int active_quad_index_given, + const unsigned int face_type) { - const auto &dof_info = data.get_dof_info(dof_no); - const auto &mapping_storage = internal::MatrixFreeFunctions:: - MappingInfoCellsOrFaces::get( - data.get_mapping_info(), quad_no); - const unsigned int active_fe_index = + typename FEEvaluationData:: + InitializationData init_data; + + init_data.dof_info = &data.get_dof_info(dof_no); + init_data.mapping_data = + &internal::MatrixFreeFunctions:: + MappingInfoCellsOrFaces::get( + data.get_mapping_info(), quad_no); + + init_data.active_fe_index = fe_degree != numbers::invalid_unsigned_int ? - dof_info.fe_index_from_degree(first_selected_component, fe_degree) : + init_data.dof_info->fe_index_from_degree(first_selected_component, + fe_degree) : (active_fe_index_given != numbers::invalid_unsigned_int ? active_fe_index_given : 0); - const unsigned int active_quad_index = + init_data.active_quad_index = fe_degree == numbers::invalid_unsigned_int ? (active_quad_index_given != numbers::invalid_unsigned_int ? active_quad_index_given : - std::min(active_fe_index, - mapping_storage.descriptor.size() - 1)) : - mapping_storage.quad_index_from_n_q_points(n_q_points); - return std::make_tuple( - &data.get_shape_info( - dof_no, - quad_no, - dof_info.component_to_base_index[first_selected_component], - active_fe_index, - active_quad_index), - &dof_info, - &mapping_storage, - active_fe_index, - active_quad_index); + std::min(init_data.active_fe_index, + init_data.mapping_data->descriptor.size() - + 1)) : + init_data.mapping_data->quad_index_from_n_q_points(n_q_points); + + init_data.shape_info = &data.get_shape_info( + dof_no, + quad_no, + init_data.dof_info->component_to_base_index[first_selected_component], + init_data.active_fe_index, + init_data.active_quad_index); + init_data.descriptor = + &init_data.mapping_data->descriptor + [is_face ? + (init_data.active_quad_index * std::max(1, dim - 1) + + (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) : + init_data.active_quad_index]; + + return init_data; } } // namespace internal @@ -2885,10 +2892,10 @@ inline FEEvaluationBase &info, - const bool is_interior_face = true, - const unsigned int quad_no = 0, - const unsigned int face_type = 0, - const unsigned int first_selected_component = 0); + FEEvaluationData(const ShapeInfoType &info, + const bool is_interior_face = true); /** * Copy constructor. @@ -436,7 +431,8 @@ public: /** * If is_face is true, this function returns the orientation index within an - * array of orientations as returned by get_orientation_map(). + * array of orientations as stored in ShapeInfo for unknowns and quadrature + * points. * * @note This function depends on the internal representation of data, which * is not stable between releases of deal.II, and is hence mostly for @@ -445,17 +441,6 @@ public: unsigned int get_face_orientation() const; - /** - * If is_face is true, this function returns the map re-ordering from face - * quadrature points in standard orientation to a specific orientation. - * - * @note This function depends on the internal representation of data, which - * is not stable between releases of deal.II, and is hence mostly for - * internal use. - */ - const Table<2, unsigned int> & - get_orientation_map() const; - /** * Return the current index in the access to compressed indices. * @@ -571,7 +556,32 @@ public: //@} + /** + * This data structure is used for the initialization by the derived + * FEEvaluation classes. + */ + struct InitializationData + { + const ShapeInfoType * shape_info; + const DoFInfo * dof_info; + const MappingInfoStorageType *mapping_data; + unsigned int active_fe_index; + unsigned int active_quad_index; + const typename MappingInfoStorageType::QuadratureDescriptor *descriptor; + }; + protected: + /** + * Constructor filling information available from a MatrixFree class, + * collected in an internal data structure to reduce overhead of + * initialization. Used in derived classes when this class is initialized + * from a MatrixFree object. + */ + FEEvaluationData(const InitializationData &info, + const bool is_interior_face, + const unsigned int quad_no, + const unsigned int first_selected_component); + /** * Constructor that comes with reduced functionality and works similarly as * FEValues. @@ -901,41 +911,46 @@ protected: template inline FEEvaluationData::FEEvaluationData( - const std::tuple &shape_dof_info, - const bool is_interior_face, - const unsigned int quad_no, - const unsigned int face_type, - const unsigned int first_selected_component) - : data(std::get<0>(shape_dof_info)) - , dof_info(std::get<1>(shape_dof_info)) - , mapping_data(std::get<2>(shape_dof_info)) + const ShapeInfoType &shape_info, + const bool is_interior_face) + : FEEvaluationData( + InitializationData{&shape_info, nullptr, nullptr, 0, 0, nullptr}, + is_interior_face, + 0, + 0) +{} + + +template +inline FEEvaluationData::FEEvaluationData( + const InitializationData &initialization_data, + const bool is_interior_face, + const unsigned int quad_no, + const unsigned int first_selected_component) + : data(initialization_data.shape_info) + , dof_info(initialization_data.dof_info) + , mapping_data(initialization_data.mapping_data) , quad_no(quad_no) , n_fe_components(dof_info != nullptr ? dof_info->start_components.back() : 0) , first_selected_component(first_selected_component) - , active_fe_index(std::get<3>(shape_dof_info)) - , active_quad_index(std::get<4>(shape_dof_info)) - , descriptor( - &mapping_data->descriptor - [is_face ? - (active_quad_index * std::max(1, dim - 1) + - (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) : - active_quad_index]) - , n_quadrature_points(descriptor->n_q_points) + , active_fe_index(initialization_data.active_fe_index) + , active_quad_index(initialization_data.active_quad_index) + , descriptor(initialization_data.descriptor) + , n_quadrature_points(is_face ? data->n_q_points_face : data->n_q_points) , jacobian(nullptr) , J_value(nullptr) , normal_vectors(nullptr) , normal_x_jacobian(nullptr) - , quadrature_weights(descriptor->quadrature_weights.begin()) + , quadrature_weights( + descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr) +# ifdef DEBUG , dof_values_initialized(false) , values_quad_initialized(false) , gradients_quad_initialized(false) , hessians_quad_initialized(false) , values_quad_submitted(false) , gradients_quad_submitted(false) +# endif , cell(numbers::invalid_unsigned_int) , is_interior_face(is_interior_face) , dof_access_index( @@ -948,7 +963,11 @@ inline FEEvaluationData::FEEvaluationData( , face_orientation(0) , subface_index(0) , cell_type(internal::MatrixFreeFunctions::general) -{} +{ + // TODO: This currently fails for simplex/matrix_free_face_integral_01 + // if (descriptor != nullptr) + // AssertDimension(n_quadrature_points, descriptor->n_q_points); +} @@ -1009,12 +1028,14 @@ FEEvaluationData::operator=(const FEEvaluationData &other) normal_x_jacobian = nullptr; quadrature_weights = other.quadrature_weights; +# ifdef DEBUG dof_values_initialized = false; values_quad_initialized = false; gradients_quad_initialized = false; hessians_quad_initialized = false; values_quad_submitted = false; gradients_quad_submitted = false; +# endif cell = numbers::invalid_unsigned_int; is_interior_face = other.is_interior_face; @@ -1132,7 +1153,7 @@ FEEvaluationData::JxW(const unsigned int q_point) const "update_values|update_gradients")); if (cell_type <= internal::MatrixFreeFunctions::affine) { - Assert(quadrature_weights != nullptr, ExcInternalError()); + Assert(quadrature_weights != nullptr, ExcNotInitialized()); return J_value[0] * quadrature_weights[q_point]; } else @@ -1397,16 +1418,6 @@ FEEvaluationData::get_face_orientation() const -template -inline const Table<2, unsigned int> & -FEEvaluationData::get_orientation_map() const -{ - Assert(descriptor != nullptr, ExcNotInitialized()); - return descriptor->face_orientations; -} - - - template inline internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex FEEvaluationData::get_dof_access_index() const diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 54fc917747..6a4c0f186b 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -126,18 +126,17 @@ namespace internal if (flag == false) { - cell_data[my_q].descriptor[hpq].initialize(quad[my_q][hpq], - update_default); + cell_data[my_q].descriptor[hpq].initialize(quad[my_q][hpq]); const auto quad_face = get_unique_face_quadratures(quad[my_q][hpq]); if (quad_face.first.size() > 0) // line, quad { face_data[my_q].descriptor[hpq * scale].initialize( - quad_face.first, update_default); + quad_face.first); face_data_by_cells[my_q] .descriptor[hpq * scale] - .initialize(quad_face.first, update_default); + .initialize(quad_face.first); } if (quad_face.second.size() > 0) // triangle @@ -145,10 +144,10 @@ namespace internal AssertDimension(dim, 3); face_data[my_q] .descriptor[hpq * scale + (is_mixed_mesh ? 1 : 0)] - .initialize(quad_face.second, update_default); + .initialize(quad_face.second); face_data_by_cells[my_q] .descriptor[hpq * scale + (is_mixed_mesh ? 1 : 0)] - .initialize(quad_face.second, update_default); + .initialize(quad_face.second); } const auto face_quadrature_collection = @@ -163,15 +162,14 @@ namespace internal else { cell_data[my_q].descriptor[hpq].initialize( - quad[my_q][hpq].get_tensor_basis()[0], update_default); + quad[my_q][hpq].get_tensor_basis()[0]); const auto quad_face = quad[my_q][hpq].get_tensor_basis()[0]; - face_data[my_q].descriptor[hpq * scale].initialize( - quad_face, this->update_flags_boundary_faces); + face_data[my_q].descriptor[hpq * scale].initialize(quad_face); 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); + quad_face); reference_cell_types[my_q][hpq] = dealii::ReferenceCells::get_hypercube(); } @@ -1131,10 +1129,7 @@ namespace internal shape_info.dofs_per_component_on_cell; constexpr unsigned int hess_dim = dim * (dim + 1) / 2; - MappingInfoStorage temp_data; - temp_data.copy_descriptor(my_data); - FEEvaluationData eval( - std::make_tuple(&shape_info, nullptr, &temp_data, 0, 0)); + FEEvaluationData eval(shape_info); AlignedVector evaluation_data; eval.set_data_pointers(&evaluation_data, dim); @@ -2040,23 +2035,23 @@ namespace internal shape_info.dofs_per_component_on_cell; constexpr unsigned int hess_dim = dim * (dim + 1) / 2; - MappingInfoStorage temp_data; - temp_data.copy_descriptor(my_data); - FEEvaluationData eval_int( - std::make_tuple(&shape_info, nullptr, &temp_data, 0, 0), true); - FEEvaluationData eval_ext( - std::make_tuple(&shape_info, nullptr, &temp_data, 0, 0), false); + FEEvaluationData eval_int(shape_info, + true); + FEEvaluationData eval_ext(shape_info, + false); - AlignedVector evaluation_data, evaluation_data_ext; + // Let both evaluators use the same array as their use will not + // overlap + AlignedVector evaluation_data; eval_int.set_data_pointers(&evaluation_data, dim); - eval_ext.set_data_pointers(&evaluation_data_ext, dim); + eval_ext.set_data_pointers(&evaluation_data, dim); for (unsigned int face = begin_face; face < end_face; ++face) for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d) { internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedDouble::size()> - face_double = {}; + face_double = {}; face_double.interior_face_no = faces[face].interior_face_no; face_double.exterior_face_no = faces[face].exterior_face_no; face_double.face_orientation = faces[face].face_orientation; diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h index 9dd239683c..98c8a3cdc1 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.h +++ b/include/deal.II/matrix_free/mapping_info_storage.h @@ -131,15 +131,13 @@ namespace internal */ template void - initialize(const Quadrature &quadrature, - const UpdateFlags update_flags_inner_faces = update_default); + initialize(const Quadrature &quadrature); /** * Set up the lengths in the various members of this struct. */ void - initialize(const Quadrature<1> &quadrature_1d, - const UpdateFlags update_flags_inner_faces = update_default); + initialize(const Quadrature<1> &quadrature_1d); /** * Returns the memory consumption in bytes. @@ -176,14 +174,6 @@ namespace internal * when it is used in a vectorized context). */ AlignedVector quadrature_weights; - - /** - * For quadrature on faces, the evaluation of basis functions is not - * in the correct order if a face is not in the standard orientation - * to a given element. This data structure is used to re-order the - * data evaluated on quadrature points to represent the correct order. - */ - dealii::Table<2, unsigned int> face_orientations; }; /** @@ -303,15 +293,6 @@ namespace internal unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const; - /** - * Get the descriptor from another instance of QuadratureDescriptor, - * clearing all other data fields. - */ - template - void - copy_descriptor( - const MappingInfoStorage &other); - /** * Helper function to determine which update flags must be set in the * internal functions to initialize all data as requested by the user. diff --git a/include/deal.II/matrix_free/mapping_info_storage.templates.h b/include/deal.II/matrix_free/mapping_info_storage.templates.h index d4c15ef09a..9a7a035000 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.templates.h +++ b/include/deal.II/matrix_free/mapping_info_storage.templates.h @@ -52,12 +52,8 @@ namespace internal template void MappingInfoStorage::QuadratureDescriptor:: - initialize(const Quadrature &quadrature, - const UpdateFlags update_flags_inner_faces) + initialize(const Quadrature &quadrature) { - Assert(structdim + 1 <= spacedim || - update_flags_inner_faces == update_default, - ExcMessage("Volume cells do not allow for setting inner faces")); this->quadrature = quadrature; n_q_points = quadrature.size(); quadrature_weights.resize(n_q_points); @@ -65,9 +61,6 @@ namespace internal quadrature_weights[i] = quadrature.weight(i); // note: quadrature_1d and tensor_quadrature_weights are not set up - - // TODO: set up face_orientations - (void)update_flags_inner_faces; } @@ -75,12 +68,8 @@ namespace internal template void MappingInfoStorage::QuadratureDescriptor:: - initialize(const Quadrature<1> &quadrature_1d, - const UpdateFlags update_flags_inner_faces) + initialize(const Quadrature<1> &quadrature_1d) { - Assert(structdim + 1 <= spacedim || - update_flags_inner_faces == update_default, - ExcMessage("Volume cells do not allow for setting inner faces")); this->quadrature_1d = quadrature_1d; quadrature = Quadrature(quadrature_1d); n_q_points = quadrature.size(); @@ -94,34 +83,6 @@ namespace internal for (unsigned int i = 0; i < quadrature_1d.size(); ++i) tensor_quadrature_weights[d][i] = quadrature_1d.weight(i); } - - // face orientation for faces in 3D - if (structdim == spacedim - 1 && spacedim == 3 && - update_flags_inner_faces != update_default) - { - const unsigned int n = quadrature_1d.size(); - face_orientations.reinit(8, n * n); - for (unsigned int j = 0, i = 0; j < n; ++j) - for (unsigned int k = 0; k < n; ++k, ++i) - { - // face_orientation=true, face_flip=false, face_rotation=false - face_orientations[0][i] = i; - // face_orientation=false, face_flip=false, face_rotation=false - face_orientations[1][i] = j + k * n; - // face_orientation=true, face_flip=true, face_rotation=false - face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; - // face_orientation=false, face_flip=true, face_rotation=false - face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; - // face_orientation=true, face_flip=false, face_rotation=true - face_orientations[4][i] = j + (n - 1 - k) * n; - // face_orientation=false, face_flip=false, face_rotation=true - face_orientations[5][i] = k + (n - 1 - j) * n; - // face_orientation=true, face_flip=true, face_rotation=true - face_orientations[6][i] = (n - 1 - j) + k * n; - // face_orientation=false, face_flip=true, face_rotation=true - face_orientations[7][i] = (n - 1 - k) + j * n; - } - } } @@ -132,8 +93,7 @@ namespace internal memory_consumption() const { std::size_t memory = sizeof(this) + quadrature.memory_consumption() + - quadrature_weights.memory_consumption() + - face_orientations.memory_consumption(); + quadrature_weights.memory_consumption(); for (int d = 0; d < structdim; ++d) memory += tensor_quadrature_weights[d].memory_consumption(); return memory; @@ -160,41 +120,6 @@ namespace internal - template - template - void - MappingInfoStorage::copy_descriptor( - const MappingInfoStorage &other) - { - clear_data_fields(); - descriptor.clear(); - descriptor.resize(other.descriptor.size()); - for (unsigned int i = 0; i < descriptor.size(); ++i) - { - descriptor[i].n_q_points = other.descriptor[i].n_q_points; - descriptor[i].quadrature_1d = other.descriptor[i].quadrature_1d; - descriptor[i].quadrature = other.descriptor[i].quadrature; - for (unsigned int d = 0; d < structdim; ++d) - { - descriptor[i].tensor_quadrature_weights[d].resize( - other.descriptor[i].tensor_quadrature_weights[d].size()); - std::copy( - other.descriptor[i].tensor_quadrature_weights[d].begin(), - other.descriptor[i].tensor_quadrature_weights[d].end(), - descriptor[i].tensor_quadrature_weights[d].begin()); - } - descriptor[i].quadrature_weights.resize( - other.descriptor[i].quadrature_weights.size()); - std::copy(other.descriptor[i].quadrature_weights.begin(), - other.descriptor[i].quadrature_weights.end(), - descriptor[i].quadrature_weights.begin()); - descriptor[i].face_orientations = - other.descriptor[i].face_orientations; - } - } - - - template UpdateFlags MappingInfoStorage::compute_update_flags( diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 241cbe517f..d7269738da 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -530,12 +530,20 @@ namespace internal dealii::Table<2, unsigned int> face_to_cell_index_hermite; /** - * For degrees on faces, the basis functions are not + * For unknowns located on faces, the basis functions are not * in the correct order if a face is not in the standard orientation * to a given element. This data structure is used to re-order the * basis functions to represent the correct order. */ - dealii::Table<2, unsigned int> face_orientations; + dealii::Table<2, unsigned int> face_orientations_dofs; + + /** + * For interpretation of values at quadrature points, the order of + * points is not correct if a face is not in the standard orientation to + * a given element. This data structure is used to re-order the + * quadrature points to represent the correct order. + */ + dealii::Table<2, unsigned int> face_orientations_quad; private: /** diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 729c7a8eb0..3d59bb92c9 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -854,40 +854,46 @@ namespace internal // (similar to MappingInfoStorage::QuadratureDescriptor::initialize) if (dim == 3) { - const unsigned int n = fe_degree + 1; - face_orientations.reinit(8, n * n); - for (unsigned int j = 0, i = 0; j < n; ++j) - for (unsigned int k = 0; k < n; ++k, ++i) - { - // face_orientation=true, face_flip=false, - // face_rotation=false - face_orientations[0][i] = i; - // face_orientation=false, face_flip=false, - // face_rotation=false - face_orientations[1][i] = j + k * n; - // face_orientation=true, face_flip=true, - // face_rotation=false - face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; - // face_orientation=false, face_flip=true, - // face_rotation=false - face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; - // face_orientation=true, face_flip=false, - // face_rotation=true - face_orientations[4][i] = j + (n - 1 - k) * n; - // face_orientation=false, face_flip=false, - // face_rotation=true - face_orientations[5][i] = k + (n - 1 - j) * n; - // face_orientation=true, face_flip=true, - // face_rotation=true - face_orientations[6][i] = (n - 1 - j) + k * n; - // face_orientation=false, face_flip=true, - // face_rotation=true - face_orientations[7][i] = (n - 1 - k) + j * n; - } + auto compute_orientations = + [](const unsigned int n, + Table<2, unsigned int> &face_orientations) { + face_orientations.reinit(8, n * n); + for (unsigned int j = 0, i = 0; j < n; ++j) + for (unsigned int k = 0; k < n; ++k, ++i) + { + // face_orientation=true, face_flip=false, + // face_rotation=false + face_orientations[0][i] = i; + // face_orientation=false, face_flip=false, + // face_rotation=false + face_orientations[1][i] = j + k * n; + // face_orientation=true, face_flip=true, + // face_rotation=false + face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; + // face_orientation=false, face_flip=true, + // face_rotation=false + face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; + // face_orientation=true, face_flip=false, + // face_rotation=true + face_orientations[4][i] = j + (n - 1 - k) * n; + // face_orientation=false, face_flip=false, + // face_rotation=true + face_orientations[5][i] = k + (n - 1 - j) * n; + // face_orientation=true, face_flip=true, + // face_rotation=true + face_orientations[6][i] = (n - 1 - j) + k * n; + // face_orientation=false, face_flip=true, + // face_rotation=true + face_orientations[7][i] = (n - 1 - k) + j * n; + } + }; + compute_orientations(fe_degree + 1, face_orientations_dofs); + compute_orientations(n_q_points_1d, face_orientations_quad); } else { - face_orientations.reinit(1, 1); + face_orientations_dofs.reinit(1, 1); + face_orientations_quad.reinit(1, 1); } }