From: Maximilian Bergbauer Date: Wed, 5 Apr 2023 15:26:44 +0000 (+0200) Subject: Optimize computations in MappingQ X-Git-Tag: v9.5.0-rc1~316^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15035%2Fhead;p=dealii.git Optimize computations in MappingQ --- diff --git a/include/deal.II/base/qprojector.h b/include/deal.II/base/qprojector.h index 381f8aa0fa..4e4749752b 100644 --- a/include/deal.II/base/qprojector.h +++ b/include/deal.II/base/qprojector.h @@ -109,6 +109,20 @@ public: const SubQuadrature &quadrature, const unsigned int face_no); + /** + * Compute the cell quadrature formula corresponding to using + * quadrature on face face_no taking into account the + * orientation of the face. For further details, see the general doc for this + * class. + */ + static Quadrature + project_to_oriented_face(const ReferenceCell &reference_cell, + const SubQuadrature &quadrature, + const unsigned int face_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation); + /** * Compute the quadrature points on the cell if the given quadrature formula * is used on face face_no, subface number subface_no @@ -148,6 +162,25 @@ public: const RefinementCase &ref_case = RefinementCase::isotropic_refinement); + /** + * Compute the cell quadrature formula corresponding to using + * quadrature on subface subface_no of face + * face_no with SubfaceCase ref_case. The last + * argument is only used in 3d. + * + * @note Only the points are transformed. The quadrature weights are the + * same as those of the original rule. + */ + static Quadrature + project_to_oriented_subface(const ReferenceCell & reference_cell, + const SubQuadrature & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const internal::SubfaceCase ref_case); + /** * Take a collection of face quadrature formulas and generate a cell * quadrature formula from it where the quadrature points of the given diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 98366d47d4..28f34de1bf 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -255,6 +255,37 @@ public: internal::FEValuesImplementation::MappingRelatedData &output_data) const; + /** + * As opposed to the fill_fe_face_values() + * function that relies on pre-computed information of InternalDataBase, this + * function chooses the flexible evaluation path on the cell and points + * passed in to the current function. + * + * @param[in] cell The cell where to evaluate the mapping. + * + * @param[in] face_number The face number where to evaluate the mapping. + * + * @param[in] face_quadrature The quadrature points where the + * transformation (Jacobians, positions) should be computed. + * + * @param[in] internal_data A reference to an object previously created + * that may be used to store information the mapping can compute once on the + * reference cell. See the documentation of the Mapping::InternalDataBase + * class for an extensive description of the purpose of these objects. + * + * @param[out] output_data A struct containing the evaluated quantities such + * as the Jacobian resulting from application of the mapping on the given + * cell with its underlying manifolds. + */ + void + fill_mapping_data_for_face_quadrature( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_number, + const Quadrature & face_quadrature, + const typename Mapping::InternalDataBase & internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const; + /** * @name Interface with FEValues and friends * @{ @@ -391,6 +422,13 @@ public: virtual std::size_t memory_consumption() const override; + /** + * Location of quadrature points of faces or subfaces in 3d with all + * possible orientations. Can be accessed with the correct offset provided + * via QProjector::DataSetDescriptor. Not needed/used for cells. + */ + AlignedVector> quadrature_points; + /** * Values of shape functions. Access by function @p shape. * diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index ddd4c145b2..e144e33096 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -1116,6 +1116,8 @@ namespace internal const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, std::vector> & quadrature_points, + std::vector> &jacobians, + std::vector> &inverse_jacobians, std::vector> &jacobian_grads) { const UpdateFlags update_flags = data.update_each; @@ -1246,6 +1248,26 @@ namespace internal data.volume_elements[point] = data.contravariant[point].determinant(); + // copy values from InternalData to vector given by reference + if (update_flags & update_jacobians) + { + const unsigned int n_q_points = data.contravariant.size(); + AssertDimension(jacobians.size(), n_q_points); + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point = 0; point < n_q_points; ++point) + jacobians[point] = data.contravariant[point]; + } + + // copy values from InternalData to vector given by reference + if (update_flags & update_inverse_jacobians) + { + const unsigned int n_q_points = data.contravariant.size(); + AssertDimension(inverse_jacobians.size(), n_q_points); + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point = 0; point < n_q_points; ++point) + inverse_jacobians[point] = data.covariant[point].transpose(); + } + if (evaluation_flag & EvaluationFlags::hessians) { constexpr int desymmetrize_3d[6][2] = { @@ -1286,112 +1308,114 @@ namespace internal } - /** - * Compute the locations of quadrature points on the object described by - * the first argument (and the cell for which the mapping support points - * have already been set), but only if the update_flags of the @p data - * argument indicate so. - */ + template inline void - maybe_compute_q_points( - const typename QProjector::DataSetDescriptor data_set, + maybe_update_q_points_Jacobians_generic( + const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, - std::vector> &quadrature_points) + const ArrayView> & unit_points, + const std::vector> &polynomials_1d, + const unsigned int polynomial_degree, + const std::vector &renumber_lexicographic_to_hierarchic, + std::vector> & quadrature_points, + std::vector> &jacobians, + std::vector> &inverse_jacobians) { - const UpdateFlags update_flags = data.update_each; + const UpdateFlags update_flags = data.update_each; + const std::vector> &support_points = + data.mapping_support_points; - if (update_flags & update_quadrature_points) - for (unsigned int point = 0; point < quadrature_points.size(); ++point) + const unsigned int n_points = unit_points.size(); + const unsigned int n_lanes = VectorizedArray::size(); + + // Use the more heavy VectorizedArray code path if there is more than + // one point left to compute + for (unsigned int i = 0; i < n_points; i += n_lanes) + if (n_points - i > 1) { - const double * shape = &data.shape(point + data_set, 0); - Point result = - (shape[0] * data.mapping_support_points[0]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - result[i] += shape[k] * data.mapping_support_points[k][i]; - quadrature_points[point] = result; - } - } + Point> p_vec; + for (unsigned int j = 0; j < n_lanes; ++j) + if (i + j < n_points) + for (unsigned int d = 0; d < dim; ++d) + p_vec[d][j] = unit_points[i + j][d]; + else + for (unsigned int d = 0; d < dim; ++d) + p_vec[d][j] = unit_points[i][d]; + const auto result = + internal::evaluate_tensor_product_value_and_gradient( + polynomials_1d, + support_points, + p_vec, + polynomial_degree == 1, + renumber_lexicographic_to_hierarchic); + if (update_flags & update_quadrature_points) + for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) + for (unsigned int d = 0; d < spacedim; ++d) + quadrature_points[i + j][d] = result.first[d][j]; - /** - * Update the co- and contravariant matrices as well as their determinant, - * for the cell - * described stored in the data object, but only if the update_flags of the @p data - * argument indicate so. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - inline void - maybe_update_Jacobians( - const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQ::InternalData &data) - { - const UpdateFlags update_flags = data.update_each; + if (cell_similarity == CellSimilarity::translation) + continue; - if (update_flags & update_contravariant_transformation) - // if the current cell is just a - // translation of the previous one, no - // need to recompute jacobians... - if (cell_similarity != CellSimilarity::translation) + if (update_flags & update_contravariant_transformation) + for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) + for (unsigned int d = 0; d < spacedim; ++d) + for (unsigned int e = 0; e < dim; ++e) + data.contravariant[i + j][d][e] = result.second[e][d][j]; + + if (update_flags & update_volume_elements) + for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) + data.volume_elements[i + j] = + data.contravariant[i + j].determinant(); + + if (update_flags & update_jacobians) + for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) + jacobians[i + j] = data.contravariant[i + j]; + + if (update_flags & update_covariant_transformation) + for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) + data.covariant[i + j] = + data.contravariant[i + j].covariant_form(); + + if (update_flags & update_inverse_jacobians) + for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) + inverse_jacobians[i + j] = data.covariant[i + j].transpose(); + } + else { - const unsigned int n_q_points = data.contravariant.size(); + const auto result = + internal::evaluate_tensor_product_value_and_gradient( + polynomials_1d, + support_points, + unit_points[i], + polynomial_degree == 1, + renumber_lexicographic_to_hierarchic); - std::fill(data.contravariant.begin(), - data.contravariant.end(), - DerivativeForm<1, dim, spacedim>()); + if (update_flags & update_quadrature_points) + quadrature_points[i] = result.first; - Assert(data.n_shape_functions > 0, ExcInternalError()); + if (cell_similarity == CellSimilarity::translation) + continue; - for (unsigned int point = 0; point < n_q_points; ++point) + if (update_flags & update_contravariant_transformation) { - double result[spacedim][dim]; - - // peel away part of sum to avoid zeroing the - // entries and adding for the first time - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - result[i][j] = data.derivative(point + data_set, 0)[j] * - data.mapping_support_points[0][i]; - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - result[i][j] += data.derivative(point + data_set, k)[j] * - data.mapping_support_points[k][i]; - - // write result into contravariant data. for - // j=dim in the case dim jac_transposed = result.second; + data.contravariant[i] = jac_transposed.transpose(); } - } - if (update_flags & update_covariant_transformation) - if (cell_similarity != CellSimilarity::translation) - { - const unsigned int n_q_points = data.contravariant.size(); - for (unsigned int point = 0; point < n_q_points; ++point) - { - data.covariant[point] = - (data.contravariant[point]).covariant_form(); - } - } + if (update_flags & update_volume_elements) + data.volume_elements[i] = data.contravariant[i].determinant(); - if (update_flags & update_volume_elements) - if (cell_similarity != CellSimilarity::translation) - { - const unsigned int n_q_points = data.contravariant.size(); - for (unsigned int point = 0; point < n_q_points; ++point) - data.volume_elements[point] = - data.contravariant[point].determinant(); + if (update_flags & update_jacobians) + jacobians[i] = data.contravariant[i]; + + if (update_flags & update_covariant_transformation) + data.covariant[i] = data.contravariant[i].covariant_form(); + + if (update_flags & update_inverse_jacobians) + inverse_jacobians[i] = data.covariant[i].transpose(); } } @@ -2012,19 +2036,11 @@ namespace internal output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / output_data.boundary_forms[i].norm()); - - if (update_flags & update_jacobians) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.jacobians[point] = data.contravariant[point]; - - if (update_flags & update_inverse_jacobians) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.inverse_jacobians[point] = - data.covariant[point].transpose(); } } + /** * Do the work of MappingQ::fill_fe_face_values() and * MappingQ::fill_fe_subface_values() in a generic way, @@ -2041,6 +2057,9 @@ namespace internal const typename QProjector::DataSetDescriptor data_set, const Quadrature & quadrature, const typename dealii::MappingQ::InternalData &data, + const std::vector> &polynomials_1d, + const unsigned int polynomial_degree, + const std::vector &renumber_lexicographic_to_hierarchic, internal::FEValuesImplementation::MappingRelatedData &output_data) { @@ -2050,16 +2069,25 @@ namespace internal CellSimilarity::none, data, output_data.quadrature_points, + output_data.jacobians, + output_data.inverse_jacobians, output_data.jacobian_grads); } else { - maybe_compute_q_points(data_set, - data, - output_data.quadrature_points); - maybe_update_Jacobians(CellSimilarity::none, - data_set, - data); + internal::MappingQImplementation:: + maybe_update_q_points_Jacobians_generic( + CellSimilarity::none, + data, + make_array_view(&data.quadrature_points[data_set], + &data.quadrature_points[data_set] + + quadrature.size()), + polynomials_1d, + polynomial_degree, + renumber_lexicographic_to_hierarchic, + output_data.quadrature_points, + output_data.jacobians, + output_data.inverse_jacobians); maybe_update_jacobian_grads( CellSimilarity::none, data_set, data, output_data.jacobian_grads); } diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index 5dc9d8f665..9dc630e15f 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -212,10 +212,11 @@ namespace NonMatching * class. */ void - compute_mapping_data_for_generic_points( + compute_mapping_data_for_quadrature( const typename Triangulation::cell_iterator &cell, - const ArrayView> & unit_points, - MappingData & mapping_data); + CellSimilarity::Similarity &cell_similarity, + const Quadrature & quadrature, + MappingData & mapping_data); /** * Compute the mapping related data for the given @p mapping, @@ -250,6 +251,12 @@ namespace NonMatching */ std::vector unit_points_index; + /** + * A pointer to the internal data of the underlying mapping. + */ + std::unique_ptr::InternalDataBase> + internal_mapping_data; + /** * A pointer to the underlying mapping. */ @@ -318,6 +325,16 @@ namespace NonMatching // always save quadrature points for now update_flags_mapping |= update_quadrature_points; + + // construct internal_mapping_data for MappingQ to be able to reuse it in + // reinit() calls to avoid memory allocations + if (const MappingQ *mapping_q = + dynamic_cast *>(&mapping)) + { + internal_mapping_data = + std::make_unique::InternalData>( + mapping_q->get_degree()); + } } @@ -328,7 +345,7 @@ namespace NonMatching const typename Triangulation::cell_iterator &cell, const std::vector> & unit_points_in) { - reinit(cell, make_array_view(unit_points_in.begin(), unit_points_in.end())); + reinit(cell, Quadrature(unit_points_in)); } @@ -339,16 +356,9 @@ namespace NonMatching const typename Triangulation::cell_iterator &cell, const ArrayView> & unit_points_in) { - unit_points = - std::vector>(unit_points_in.begin(), unit_points_in.end()); - - mapping_data.resize(1); - compute_mapping_data_for_generic_points(cell, - unit_points_in, - mapping_data[0]); - - state = State::single_cell; - is_reinitialized(); + reinit(cell, + std::vector>(unit_points_in.begin(), + unit_points_in.end())); } @@ -359,16 +369,18 @@ namespace NonMatching const typename Triangulation::cell_iterator &cell, const Quadrature & quadrature) { - const auto &points = quadrature.get_points(); - const auto &weights = quadrature.get_weights(); + unit_points = quadrature.get_points(); - reinit(cell, points); + mapping_data.resize(1); + CellSimilarity::Similarity cell_similarity = + CellSimilarity::Similarity::none; + compute_mapping_data_for_quadrature(cell, + cell_similarity, + quadrature, + mapping_data[0]); - if (update_flags_mapping & update_JxW_values) - for (unsigned int q = 0; q < points.size(); ++q) - mapping_data[0].JxW_values[q] = - determinant(Tensor<2, dim>(mapping_data[0].jacobians[q])) * - weights[q]; + state = State::single_cell; + is_reinitialized(); } @@ -380,11 +392,34 @@ namespace NonMatching const IteratorRange & cell_iterator_range, const std::vector>> &unit_points_vector, const unsigned int n_unfiltered_cells) + { + const unsigned int n_cells = unit_points_vector.size(); + AssertDimension(n_cells, + std::distance(cell_iterator_range.begin(), + cell_iterator_range.end())); + + std::vector> quadrature_vector(n_cells); + for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index) + quadrature_vector[cell_index] = + Quadrature(quadrature_vector[cell_index].get_points()); + + reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells); + } + + + + template + template + void + MappingInfo::reinit_cells( + const IteratorRange & cell_iterator_range, + const std::vector> &quadrature_vector, + const unsigned int n_unfiltered_cells) { do_cell_index_compression = n_unfiltered_cells != numbers::invalid_unsigned_int; - const unsigned int n_cells = unit_points_vector.size(); + const unsigned int n_cells = quadrature_vector.size(); AssertDimension(n_cells, std::distance(cell_iterator_range.begin(), cell_iterator_range.end())); @@ -392,9 +427,8 @@ namespace NonMatching // fill unit points index offset vector unit_points_index.reserve(n_cells + 1); unit_points_index.push_back(0); - for (const auto &unit_points : unit_points_vector) - unit_points_index.push_back(unit_points_index.back() + - unit_points.size()); + for (const auto &quadrature : quadrature_vector) + unit_points_index.push_back(unit_points_index.back() + quadrature.size()); const unsigned int n_unit_points = unit_points_index.back(); @@ -404,19 +438,23 @@ namespace NonMatching if (do_cell_index_compression) cell_index_to_compressed_cell_index.resize(n_unfiltered_cells, numbers::invalid_unsigned_int); + CellSimilarity::Similarity cell_similarity = + CellSimilarity::Similarity::none; unsigned int cell_index = 0; for (const auto &cell : cell_iterator_range) { auto it = unit_points.begin() + unit_points_index[cell_index]; - for (const auto &unit_point : unit_points_vector[cell_index]) + for (const auto &unit_point : + quadrature_vector[cell_index].get_points()) { *it = unit_point; ++it; } - compute_mapping_data_for_generic_points(cell, - unit_points_vector[cell_index], - mapping_data[cell_index]); + compute_mapping_data_for_quadrature(cell, + cell_similarity, + quadrature_vector[cell_index], + mapping_data[cell_index]); if (do_cell_index_compression) cell_index_to_compressed_cell_index[cell->active_cell_index()] = @@ -431,41 +469,6 @@ namespace NonMatching - template - template - void - MappingInfo::reinit_cells( - const IteratorRange & cell_iterator_range, - const std::vector> &quadrature_vector, - const unsigned int n_unfiltered_cells) - { - const unsigned int n_cells = quadrature_vector.size(); - AssertDimension(n_cells, - std::distance(cell_iterator_range.begin(), - cell_iterator_range.end())); - - std::vector>> unit_points_vector(n_cells); - for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index) - unit_points_vector[cell_index] = std::vector>( - quadrature_vector[cell_index].get_points().begin(), - quadrature_vector[cell_index].get_points().end()); - - reinit_cells(cell_iterator_range, unit_points_vector, n_unfiltered_cells); - - if (update_flags_mapping & update_JxW_values) - for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index) - { - const auto &weights = quadrature_vector[cell_index].get_weights(); - for (unsigned int q = 0; q < weights.size(); ++q) - mapping_data[cell_index].JxW_values[q] = - determinant( - Tensor<2, dim>(mapping_data[cell_index].jacobians[q])) * - weights[q]; - } - } - - - template template void @@ -809,47 +812,35 @@ namespace NonMatching template void - MappingInfo::compute_mapping_data_for_generic_points( + MappingInfo::compute_mapping_data_for_quadrature( const typename Triangulation::cell_iterator &cell, - const ArrayView> & unit_points, + CellSimilarity::Similarity & cell_similarity, + const Quadrature & quadrature, MappingData & mapping_data) { + update_flags_mapping |= + mapping->requires_update_flags(update_flags_mapping); + + mapping_data.initialize(quadrature.size(), update_flags_mapping); + + // reuse internal_mapping_data for MappingQ to avoid memory allocations if (const MappingQ *mapping_q = dynamic_cast *>(&(*mapping))) { - mapping_q->fill_mapping_data_for_generic_points(cell, - unit_points, - update_flags_mapping, - mapping_data); - } - else if (const MappingCartesian *mapping_cartesian = - dynamic_cast *>( - &(*mapping))) - { - mapping_cartesian->fill_mapping_data_for_generic_points( - cell, unit_points, update_flags_mapping, mapping_data); + (void)mapping_q; + auto &data = + dynamic_cast::InternalData &>( + *internal_mapping_data); + data.initialize(update_flags_mapping, quadrature, quadrature.size()); } else { - FE_DGQ dummy_fe(1); - dealii::FEValues fe_values( - *mapping, - dummy_fe, - Quadrature( - std::vector>(unit_points.begin(), unit_points.end())), - update_flags_mapping); - fe_values.reinit(cell); - mapping_data.initialize(unit_points.size(), update_flags_mapping); - if (update_flags_mapping & update_jacobians) - for (unsigned int q = 0; q < unit_points.size(); ++q) - mapping_data.jacobians[q] = fe_values.jacobian(q); - if (update_flags_mapping & update_inverse_jacobians) - for (unsigned int q = 0; q < unit_points.size(); ++q) - mapping_data.inverse_jacobians[q] = fe_values.inverse_jacobian(q); - if (update_flags_mapping & update_quadrature_points) - for (unsigned int q = 0; q < unit_points.size(); ++q) - mapping_data.quadrature_points[q] = fe_values.quadrature_point(q); + internal_mapping_data = + mapping->get_data(update_flags_mapping, quadrature); } + + cell_similarity = mapping->fill_fe_values( + cell, cell_similarity, quadrature, *internal_mapping_data, mapping_data); } @@ -865,11 +856,23 @@ namespace NonMatching update_flags_mapping |= mapping->requires_update_flags(update_flags_mapping); - mapping_data.initialize(quadrature.get_points().size(), - update_flags_mapping); + mapping_data.initialize(quadrature.size(), update_flags_mapping); - auto internal_mapping_data = - mapping->get_data(update_flags_mapping, quadrature); + // reuse internal_mapping_data for MappingQ to avoid memory allocations + if (const MappingQ *mapping_q = + dynamic_cast *>(&(*mapping))) + { + (void)mapping_q; + auto &data = + dynamic_cast::InternalData &>( + *internal_mapping_data); + data.initialize(update_flags_mapping, quadrature, quadrature.size()); + } + else + { + internal_mapping_data = + mapping->get_data(update_flags_mapping, quadrature); + } mapping->fill_fe_immersed_surface_values(cell, quadrature, @@ -890,18 +893,40 @@ namespace NonMatching update_flags_mapping |= mapping->requires_update_flags(update_flags_mapping); - mapping_data.initialize(quadrature.get_points().size(), - update_flags_mapping); - - auto internal_mapping_data = - mapping->get_face_data(update_flags_mapping, - hp::QCollection(quadrature)); + mapping_data.initialize(quadrature.size(), update_flags_mapping); - mapping->fill_fe_face_values(cell, - face_no, - hp::QCollection(quadrature), - *internal_mapping_data, - mapping_data); + // reuse internal_mapping_data for MappingQ to avoid memory allocations + if (const MappingQ *mapping_q = + dynamic_cast *>(&(*mapping))) + { + auto &data = + dynamic_cast::InternalData &>( + *internal_mapping_data); + data.initialize_face(update_flags_mapping, + QProjector::project_to_oriented_face( + ReferenceCells::get_hypercube(), + quadrature, + face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no)), + quadrature.size()); + + mapping_q->fill_mapping_data_for_face_quadrature( + cell, face_no, quadrature, *internal_mapping_data, mapping_data); + } + else + { + auto internal_mapping_data = + mapping->get_face_data(update_flags_mapping, + hp::QCollection(quadrature)); + + mapping->fill_fe_face_values(cell, + face_no, + hp::QCollection(quadrature), + *internal_mapping_data, + mapping_data); + } } } // namespace NonMatching diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 92600865dc..e96aed68ec 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -143,6 +143,193 @@ namespace internal q_points.push_back(cell_point); } } + + std::vector> + mutate_points_with_offset(const std::vector> &points, + const unsigned int offset) + { + switch (offset) + { + case 0: + return points; + case 1: + case 2: + case 3: + return rotate(points, offset); + case 4: + return reflect(points); + case 5: + case 6: + case 7: + return rotate(reflect(points), 8 - offset); + default: + Assert(false, ExcInternalError()); + } + return {}; + } + + Quadrature<2> + mutate_quadrature(const Quadrature<2> &quadrature, + const bool face_orientation, + const bool face_flip, + const bool face_rotation) + { + static const unsigned int offset[2][2][2] = { + {{4, 5}, // face_orientation=false; face_flip=false; + // face_rotation=false and true + {6, 7}}, // face_orientation=false; face_flip=true; + // face_rotation=false and true + {{0, 1}, // face_orientation=true; face_flip=false; + // face_rotation=false and true + {2, 3}}}; // face_orientation=true; face_flip=true; + // face_rotation=false and true + + return Quadrature<2>( + mutate_points_with_offset( + quadrature.get_points(), + offset[face_orientation][face_flip][face_rotation]), + quadrature.get_weights()); + } + + std::pair> + select_subface_no_and_refinement_case( + const unsigned int subface_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const internal::SubfaceCase<3> ref_case) + { + constexpr int dim = 3; + // for each subface of a given FaceRefineCase + // there is a corresponding equivalent + // subface number of one of the "standard" + // RefineCases (cut_x, cut_y, cut_xy). Map + // the given values to those equivalent + // ones. + + // first, define an invalid number + static const unsigned int e = numbers::invalid_unsigned_int; + + static const RefinementCase + equivalent_refine_case[internal::SubfaceCase::case_isotropic + 1] + [GeometryInfo<3>::max_children_per_face] = { + // case_none. there should be only + // invalid values here. However, as + // this function is also called (in + // tests) for cells which have no + // refined faces, use isotropic + // refinement instead + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_x + {RefinementCase::cut_x, + RefinementCase::cut_x, + RefinementCase::no_refinement, + RefinementCase::no_refinement}, + // case_x1y + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_x, + RefinementCase::no_refinement}, + // case_x2y + {RefinementCase::cut_x, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::no_refinement}, + // case_x1y2y + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_y + {RefinementCase::cut_y, + RefinementCase::cut_y, + RefinementCase::no_refinement, + RefinementCase::no_refinement}, + // case_y1x + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_y, + RefinementCase::no_refinement}, + // case_y2x + {RefinementCase::cut_y, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::no_refinement}, + // case_y1x2x + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_xy (case_isotropic) + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}}; + + static const unsigned int + equivalent_subface_number[internal::SubfaceCase::case_isotropic + + 1][GeometryInfo<3>::max_children_per_face] = + {// case_none, see above + {0, 1, 2, 3}, + // case_x + {0, 1, e, e}, + // case_x1y + {0, 2, 1, e}, + // case_x2y + {0, 1, 3, e}, + // case_x1y2y + {0, 2, 1, 3}, + // case_y + {0, 1, e, e}, + // case_y1x + {0, 1, 1, e}, + // case_y2x + {0, 2, 3, e}, + // case_y1x2x + {0, 1, 2, 3}, + // case_xy (case_isotropic) + {0, 1, 2, 3}}; + + // If face-orientation or face_rotation are + // non-standard, cut_x and cut_y have to be + // exchanged. + static const RefinementCase ref_case_permutation[4] = { + RefinementCase::no_refinement, + RefinementCase::cut_y, + RefinementCase::cut_x, + RefinementCase::cut_xy}; + + // set a corresponding (equivalent) + // RefineCase and subface number + const RefinementCase equ_ref_case = + equivalent_refine_case[ref_case][subface_no]; + const unsigned int equ_subface_no = + equivalent_subface_number[ref_case][subface_no]; + // make sure, that we got a valid subface and RefineCase + Assert(equ_ref_case != RefinementCase::no_refinement, + ExcInternalError()); + Assert(equ_subface_no != e, ExcInternalError()); + // now, finally respect non-standard faces + const RefinementCase final_ref_case = + (face_orientation == face_rotation ? + ref_case_permutation[equ_ref_case] : + equ_ref_case); + + const unsigned int final_subface_no = + GeometryInfo::child_cell_on_face(RefinementCase( + final_ref_case), + 4, + equ_subface_no, + face_orientation, + face_flip, + face_rotation, + equ_ref_case); + + return std::make_pair(final_subface_no, final_ref_case); + } } // namespace } // namespace QProjector } // namespace internal @@ -250,6 +437,38 @@ QProjector<3>::project_to_face(const ReferenceCell & reference_cell, } +template +Quadrature +QProjector::project_to_oriented_face(const ReferenceCell &reference_cell, + const Quadrature &quadrature, + const unsigned int face_no, + const bool, + const bool, + const bool) +{ + return QProjector::project_to_face(reference_cell, quadrature, face_no); +} + + + +template <> +Quadrature<3> +QProjector<3>::project_to_oriented_face(const ReferenceCell &reference_cell, + const Quadrature<2> &quadrature, + const unsigned int face_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation) +{ + Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented()); + + const Quadrature<2> mutation = internal::QProjector::mutate_quadrature( + quadrature, face_orientation, face_flip, face_rotation); + + return QProjector<3>::project_to_face(reference_cell, mutation, face_no); +} + + template <> void @@ -440,6 +659,60 @@ QProjector<3>::project_to_subface(const ReferenceCell & reference_cell, +template +Quadrature +QProjector::project_to_oriented_subface( + const ReferenceCell & reference_cell, + const Quadrature &quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const bool, + const bool, + const bool, + const internal::SubfaceCase) +{ + return QProjector::project_to_subface( + reference_cell, + quadrature, + face_no, + subface_no, + RefinementCase::isotropic_refinement); +} + + + +template <> +Quadrature<3> +QProjector<3>::project_to_oriented_subface( + const ReferenceCell & reference_cell, + const Quadrature<2> & quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const internal::SubfaceCase<3> ref_case) +{ + Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented()); + + const Quadrature<2> mutation = internal::QProjector::mutate_quadrature( + quadrature, face_orientation, face_flip, face_rotation); + + const std::pair> + final_subface_no_and_ref_case = + internal::QProjector::select_subface_no_and_refinement_case( + subface_no, face_orientation, face_flip, face_rotation, ref_case); + + return QProjector<3>::project_to_subface( + reference_cell, + mutation, + face_no, + final_subface_no_and_ref_case.first, + final_subface_no_and_ref_case.second); +} + + + template <> Quadrature<1> QProjector<1>::project_to_all_faces(const ReferenceCell & reference_cell, @@ -751,29 +1024,10 @@ QProjector<3>::project_to_all_faces(const ReferenceCell & reference_cell, std::vector weights; weights.reserve(n_points_total); - Table<2, std::vector>> mutations(quadrature.size(), 8); - - // do the following for all possible mutations of a face (mutation==0 - // corresponds to a face with standard orientation, no flip and no - // rotation) - for (unsigned int face = 0; face < quadrature.size(); ++face) - { - const std::vector> &quad = quadrature[face].get_points(); - mutations(face, 0) = quad; - mutations(face, 1) = internal::QProjector::rotate(quad, 1); - mutations(face, 2) = internal::QProjector::rotate(quad, 2); - mutations(face, 3) = internal::QProjector::rotate(quad, 3); - mutations(face, 4) = internal::QProjector::reflect(quad); - mutations(face, 5) = - internal::QProjector::rotate(mutations(face, 4), 3); - mutations(face, 6) = - internal::QProjector::rotate(mutations(face, 4), 2); - mutations(face, 7) = - internal::QProjector::rotate(mutations(face, 4), 1); - } - - for (unsigned int i = 0; i < 8; ++i) + for (unsigned int offset = 0; offset < 8; ++offset) { + const auto mutation = internal::QProjector::mutate_points_with_offset( + quadrature[0].get_points(), offset); // project to each face and append results for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) @@ -781,7 +1035,11 @@ QProjector<3>::project_to_all_faces(const ReferenceCell & reference_cell, const unsigned int q_index = quadrature.size() == 1 ? 0 : face; internal::QProjector::project_to_hex_face_and_append( - mutations(q_index, i), face, q_points); + q_index > 0 ? internal::QProjector::mutate_points_with_offset( + quadrature[face].get_points(), offset) : + mutation, + face, + q_points); std::copy(quadrature[q_index].get_weights().begin(), quadrature[q_index].get_weights().end(), @@ -789,7 +1047,6 @@ QProjector<3>::project_to_all_faces(const ReferenceCell & reference_cell, } } - Assert(q_points.size() == n_points_total, ExcInternalError()); Assert(weights.size() == n_points_total, ExcInternalError()); @@ -908,18 +1165,7 @@ QProjector<3>::project_to_all_subfaces(const ReferenceCell &reference_cell, Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented()); - const unsigned int dim = 3; - const std::vector> &quad = quadrature.get_points(); - std::vector> q_reflected = internal::QProjector::reflect(quad); - std::array>, 8> mutations{ - {quad, - internal::QProjector::rotate(quad, 1), - internal::QProjector::rotate(quad, 2), - internal::QProjector::rotate(quad, 3), - q_reflected, - internal::QProjector::rotate(q_reflected, 3), - internal::QProjector::rotate(q_reflected, 2), - internal::QProjector::rotate(q_reflected, 1)}}; + const unsigned int dim = 3; const unsigned int n_points = quadrature.size(), n_faces = GeometryInfo::faces_per_cell, @@ -934,8 +1180,12 @@ QProjector<3>::project_to_all_subfaces(const ReferenceCell &reference_cell, // do the following for all possible mutations of a face (mutation==0 // corresponds to a face with standard orientation, no flip and no rotation) - for (const auto &mutation : mutations) + for (unsigned int offset = 0; offset < 8; ++offset) { + const auto mutation = + internal::QProjector::mutate_points_with_offset(quadrature.get_points(), + offset); + // project to each face and copy results for (unsigned int face = 0; face < n_faces; ++face) for (unsigned int ref_case = RefinementCase::cut_xy; @@ -1435,161 +1685,14 @@ QProjector<3>::DataSetDescriptor::subface( 0 // cut_xy }; - - // for each subface of a given FaceRefineCase - // there is a corresponding equivalent - // subface number of one of the "standard" - // RefineCases (cut_x, cut_y, cut_xy). Map - // the given values to those equivalent - // ones. - - // first, define an invalid number - static const unsigned int e = numbers::invalid_unsigned_int; - - static const RefinementCase - equivalent_refine_case[internal::SubfaceCase::case_isotropic + 1] - [GeometryInfo<3>::max_children_per_face] = { - // case_none. there should be only - // invalid values here. However, as - // this function is also called (in - // tests) for cells which have no - // refined faces, use isotropic - // refinement instead - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_x - {RefinementCase::cut_x, - RefinementCase::cut_x, - RefinementCase::no_refinement, - RefinementCase::no_refinement}, - // case_x1y - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_x, - RefinementCase::no_refinement}, - // case_x2y - {RefinementCase::cut_x, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::no_refinement}, - // case_x1y2y - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_y - {RefinementCase::cut_y, - RefinementCase::cut_y, - RefinementCase::no_refinement, - RefinementCase::no_refinement}, - // case_y1x - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_y, - RefinementCase::no_refinement}, - // case_y2x - {RefinementCase::cut_y, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::no_refinement}, - // case_y1x2x - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_xy (case_isotropic) - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}}; - - static const unsigned int - equivalent_subface_number[internal::SubfaceCase::case_isotropic + 1] - [GeometryInfo<3>::max_children_per_face] = { - // case_none, see above - {0, 1, 2, 3}, - // case_x - {0, 1, e, e}, - // case_x1y - {0, 2, 1, e}, - // case_x2y - {0, 1, 3, e}, - // case_x1y2y - {0, 2, 1, 3}, - // case_y - {0, 1, e, e}, - // case_y1x - {0, 1, 1, e}, - // case_y2x - {0, 2, 3, e}, - // case_y1x2x - {0, 1, 2, 3}, - // case_xy (case_isotropic) - {0, 1, 2, 3}}; - - // If face-orientation or face_rotation are - // non-standard, cut_x and cut_y have to be - // exchanged. - static const RefinementCase ref_case_permutation[4] = { - RefinementCase::no_refinement, - RefinementCase::cut_y, - RefinementCase::cut_x, - RefinementCase::cut_xy}; - - // set a corresponding (equivalent) - // RefineCase and subface number - const RefinementCase equ_ref_case = - equivalent_refine_case[ref_case][subface_no]; - const unsigned int equ_subface_no = - equivalent_subface_number[ref_case][subface_no]; - // make sure, that we got a valid subface and RefineCase - Assert(equ_ref_case != RefinementCase::no_refinement, - ExcInternalError()); - Assert(equ_subface_no != e, ExcInternalError()); - // now, finally respect non-standard faces - const RefinementCase final_ref_case = - (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] : - equ_ref_case); - - // what we have now is the number of - // the subface in the natural - // orientation of the *face*. what we - // need to know is the number of the - // subface concerning the standard face - // orientation as seen from the *cell*. - - // this mapping is not trivial, but we - // have done exactly this stuff in the - // child_cell_on_face function. in - // order to reduce the amount of code - // as well as to make maintaining the - // functionality easier we want to - // reuse that information. So we note - // that on the bottom face (face 4) of - // a hex cell the local x and y - // coordinates of the face and the cell - // coincide, thus also the refinement - // case of the face corresponds to the - // refinement case of the cell - // (ignoring cell refinement along the - // z direction). Using this knowledge - // we can (ab)use the - // child_cell_on_face function to do - // exactly the transformation we are in - // need of now - const unsigned int final_subface_no = - GeometryInfo::child_cell_on_face(RefinementCase(final_ref_case), - 4, - equ_subface_no, - face_orientation, - face_flip, - face_rotation, - equ_ref_case); + const std::pair> + final_subface_no_and_ref_case = + internal::QProjector::select_subface_no_and_refinement_case( + subface_no, face_orientation, face_flip, face_rotation, ref_case); return (((face_no * total_subfaces_per_face + - ref_case_offset[final_ref_case - 1] + final_subface_no) + + ref_case_offset[final_subface_no_and_ref_case.second - 1] + + final_subface_no_and_ref_case.first) + orientation_offset[face_orientation][face_flip][face_rotation]) * n_quadrature_points); } diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index d945c3a066..2961f73026 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -83,21 +83,14 @@ template void MappingQ::InternalData::initialize( const UpdateFlags update_flags, - const Quadrature &q, + const Quadrature &quadrature, const unsigned int n_original_q_points) { // store the flags in the internal data object so we can access them // in fill_fe_*_values() this->update_each = update_flags; - const unsigned int n_q_points = q.size(); - - const bool needs_higher_order_terms = - this->update_each & - (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives | - update_jacobian_pushed_forward_2nd_derivatives | - update_jacobian_3rd_derivatives | - update_jacobian_pushed_forward_3rd_derivatives); + const unsigned int n_q_points = quadrature.size(); if (this->update_each & update_covariant_transformation) covariant.resize(n_original_q_points); @@ -108,7 +101,7 @@ MappingQ::InternalData::initialize( if (this->update_each & update_volume_elements) volume_elements.resize(n_original_q_points); - tensor_product_quadrature = q.is_tensor_product(); + tensor_product_quadrature = quadrature.is_tensor_product(); // use of MatrixFree only for higher order elements and with more than one // point where tensor products do not make sense @@ -122,7 +115,7 @@ MappingQ::InternalData::initialize( if (tensor_product_quadrature) { const std::array, dim> &quad_array = - q.get_tensor_basis(); + quadrature.get_tensor_basis(); for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i) { if (quad_array[i - 1].size() != quad_array[i].size()) @@ -158,37 +151,43 @@ MappingQ::InternalData::initialize( // numbering manually (building an FE_Q is relatively // expensive due to constraints) const FE_DGQ<1> fe(polynomial_degree); - shape_info.reinit(q.get_tensor_basis()[0], fe); + shape_info.reinit(quadrature.get_tensor_basis()[0], fe); shape_info.lexicographic_numbering = FETools::lexicographic_to_hierarchic_numbering( polynomial_degree); - shape_info.n_q_points = q.size(); + shape_info.n_q_points = n_q_points; shape_info.dofs_per_component_on_cell = Utilities::pow(polynomial_degree + 1, dim); } } } + const bool needs_higher_order_terms = + this->update_each & + (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives | + update_jacobian_pushed_forward_2nd_derivatives | + update_jacobian_3rd_derivatives | + update_jacobian_pushed_forward_3rd_derivatives); + + const bool needs_higher_order_terms_generic = + !tensor_product_quadrature && + (needs_higher_order_terms || this->update_each & update_jacobian_grads); + // Only fill the big arrays on demand in case we cannot use the tensor // product quadrature code path - if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms) + if (dim == 1 || needs_higher_order_terms_generic || needs_higher_order_terms) { - // see if we need the (transformation) shape function values - // and/or gradients and resize the necessary arrays - if (this->update_each & update_quadrature_points) - shape_values.resize(n_shape_functions * n_q_points); - - if (this->update_each & - (update_covariant_transformation | - update_contravariant_transformation | update_JxW_values | - update_boundary_forms | update_normal_vectors | update_jacobians | - update_jacobian_grads | update_inverse_jacobians | - update_jacobian_pushed_forward_grads | - update_jacobian_2nd_derivatives | - update_jacobian_pushed_forward_2nd_derivatives | - update_jacobian_3rd_derivatives | - update_jacobian_pushed_forward_3rd_derivatives)) - shape_derivatives.resize(n_shape_functions * n_q_points); + // compute shapes and derivatives for codim1 (for + // do_transform_real_to_unit_cell_internal_codim1) + if (dim == (spacedim - 1)) + { + // see if we need the (transformation) shape function values + // and/or gradients and resize the necessary arrays + if (this->update_each & update_quadrature_points) + shape_values.resize(n_shape_functions * n_q_points); + if (this->update_each & update_jacobians) + shape_derivatives.resize(n_shape_functions * n_q_points); + } if (this->update_each & (update_jacobian_grads | update_jacobian_pushed_forward_grads)) @@ -203,7 +202,7 @@ MappingQ::InternalData::initialize( shape_fourth_derivatives.resize(n_shape_functions * n_q_points); // now also fill the various fields with their correct values - compute_shape_function_values(q.get_points()); + compute_shape_function_values(quadrature.get_points()); } } @@ -213,16 +212,21 @@ template void MappingQ::InternalData::initialize_face( const UpdateFlags update_flags, - const Quadrature &q, + const Quadrature &quadrature, const unsigned int n_original_q_points) { - initialize(update_flags, q, n_original_q_points); + initialize(update_flags, quadrature, n_original_q_points); + + const unsigned int n_q_points = quadrature.size(); + quadrature_points.resize(n_q_points); + for (unsigned int q = 0; q < n_q_points; ++q) + quadrature_points[q] = quadrature.get_points()[q]; if (dim > 1 && tensor_product_quadrature) { constexpr unsigned int facedim = dim - 1; const FE_DGQ<1> fe(polynomial_degree); - shape_info.reinit(q.get_tensor_basis()[0], fe); + shape_info.reinit(quadrature.get_tensor_basis()[0], fe); shape_info.lexicographic_numbering = FETools::lexicographic_to_hierarchic_numbering( polynomial_degree); @@ -985,19 +989,22 @@ MappingQ::fill_fe_values( computed_cell_similarity, data, output_data.quadrature_points, + output_data.jacobians, + output_data.inverse_jacobians, output_data.jacobian_grads); } else { - internal::MappingQImplementation::maybe_compute_q_points( - QProjector::DataSetDescriptor::cell(), - data, - output_data.quadrature_points); - - internal::MappingQImplementation::maybe_update_Jacobians( + internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic( computed_cell_similarity, - QProjector::DataSetDescriptor::cell(), - data); + data, + make_array_view(quadrature.get_points()), + polynomials_1d, + polynomial_degree, + renumber_lexicographic_to_hierarchic, + output_data.quadrature_points, + output_data.jacobians, + output_data.inverse_jacobians); internal::MappingQImplementation::maybe_update_jacobian_grads( @@ -1135,27 +1142,6 @@ MappingQ::fill_fe_values( } } - - - // copy values from InternalData to vector given by reference - if (update_flags & update_jacobians) - { - AssertDimension(output_data.jacobians.size(), n_q_points); - if (computed_cell_similarity != CellSimilarity::translation) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.jacobians[point] = data.contravariant[point]; - } - - // copy values from InternalData to vector given by reference - if (update_flags & update_inverse_jacobians) - { - AssertDimension(output_data.inverse_jacobians.size(), n_q_points); - if (computed_cell_similarity != CellSimilarity::translation) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.inverse_jacobians[point] = - data.covariant[point].transpose(); - } - return computed_cell_similarity; } @@ -1205,6 +1191,9 @@ MappingQ::fill_fe_face_values( quadrature[0].size()), quadrature[0], data, + polynomials_1d, + polynomial_degree, + renumber_lexicographic_to_hierarchic, output_data); } @@ -1255,6 +1244,9 @@ MappingQ::fill_fe_subface_values( cell->subface_case(face_no)), quadrature, data, + polynomials_1d, + polynomial_degree, + renumber_lexicographic_to_hierarchic, output_data); } @@ -1281,13 +1273,16 @@ MappingQ::fill_fe_immersed_surface_values( data.mapping_support_points = this->compute_mapping_support_points(cell); data.cell_of_current_support_points = cell; - internal::MappingQImplementation::maybe_compute_q_points( - QProjector::DataSetDescriptor::cell(), + internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic( + CellSimilarity::none, data, - output_data.quadrature_points); - - internal::MappingQImplementation::maybe_update_Jacobians( - CellSimilarity::none, QProjector::DataSetDescriptor::cell(), data); + make_array_view(quadrature.get_points()), + polynomials_1d, + polynomial_degree, + renumber_lexicographic_to_hierarchic, + output_data.quadrature_points, + output_data.jacobians, + output_data.inverse_jacobians); internal::MappingQImplementation::maybe_update_jacobian_grads( CellSimilarity::none, @@ -1372,23 +1367,6 @@ MappingQ::fill_fe_immersed_surface_values( } } } - - // copy values from InternalData to vector given by reference - if ((update_flags & update_jacobians) != 0u) - { - AssertDimension(output_data.jacobians.size(), n_q_points); - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.jacobians[point] = data.contravariant[point]; - } - - // copy values from InternalData to vector given by reference - if ((update_flags & update_inverse_jacobians) != 0u) - { - AssertDimension(output_data.inverse_jacobians.size(), n_q_points); - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.inverse_jacobians[point] = - data.covariant[point].transpose(); - } } @@ -1411,85 +1389,60 @@ MappingQ::fill_mapping_data_for_generic_points( ExcNotImplemented()); output_data.initialize(unit_points.size(), update_flags); - const std::vector> support_points = - this->compute_mapping_support_points(cell); - const unsigned int n_points = unit_points.size(); - const unsigned int n_lanes = VectorizedArray::size(); + auto internal_data = + this->get_data(update_flags, + Quadrature(std::vector>(unit_points.begin(), + unit_points.end()))); + const InternalData &data = static_cast(*internal_data); + data.mapping_support_points = this->compute_mapping_support_points(cell); - // Use the more heavy VectorizedArray code path if there is more than - // one point left to compute - for (unsigned int i = 0; i < n_points; i += n_lanes) - if (n_points - i > 1) - { - Point> p_vec; - for (unsigned int j = 0; j < n_lanes; ++j) - if (i + j < n_points) - for (unsigned int d = 0; d < dim; ++d) - p_vec[d][j] = unit_points[i + j][d]; - else - for (unsigned int d = 0; d < dim; ++d) - p_vec[d][j] = unit_points[i][d]; - - const auto result = - internal::evaluate_tensor_product_value_and_gradient( - polynomials_1d, - support_points, - p_vec, - polynomial_degree == 1, - renumber_lexicographic_to_hierarchic); - - if (update_flags & update_quadrature_points) - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - for (unsigned int d = 0; d < spacedim; ++d) - output_data.quadrature_points[i + j][d] = result.first[d][j]; + internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic( + CellSimilarity::none, + data, + unit_points, + polynomials_1d, + polynomial_degree, + renumber_lexicographic_to_hierarchic, + output_data.quadrature_points, + output_data.jacobians, + output_data.inverse_jacobians); +} - if (update_flags & update_jacobians) - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int e = 0; e < dim; ++e) - output_data.jacobians[i + j][d][e] = result.second[e][d][j]; - if (update_flags & update_inverse_jacobians) - { - DerivativeForm<1, spacedim, dim, VectorizedArray> jac( - result.second); - const DerivativeForm<1, spacedim, dim, VectorizedArray> - inv_jac = jac.covariant_form(); - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < spacedim; ++e) - output_data.inverse_jacobians[i + j][d][e] = inv_jac[d][e][j]; - } - } - else - { - const auto result = - internal::evaluate_tensor_product_value_and_gradient( - polynomials_1d, - support_points, - unit_points[i], - polynomial_degree == 1, - renumber_lexicographic_to_hierarchic); - - if (update_flags & update_quadrature_points) - output_data.quadrature_points[i] = result.first; - - if (update_flags & update_jacobians) - { - DerivativeForm<1, spacedim, dim> jac = result.second; - output_data.jacobians[i] = jac.transpose(); - } - if (update_flags & update_inverse_jacobians) - { - DerivativeForm<1, spacedim, dim> jac(result.second); - DerivativeForm<1, spacedim, dim> inv_jac = jac.covariant_form(); - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < spacedim; ++e) - output_data.inverse_jacobians[i][d][e] = inv_jac[d][e]; - } - } +template +void +MappingQ::fill_mapping_data_for_face_quadrature( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & face_quadrature, + const typename Mapping::InternalDataBase & internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + if (face_quadrature.get_points().empty()) + return; + + // ensure that the following static_cast is really correct: + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + data.mapping_support_points = this->compute_mapping_support_points(cell); + + internal::MappingQImplementation::do_fill_fe_face_values( + *this, + cell, + face_no, + numbers::invalid_unsigned_int, + QProjector::DataSetDescriptor::cell(), + face_quadrature, + data, + polynomials_1d, + polynomial_degree, + renumber_lexicographic_to_hierarchic, + output_data); }