From bb27ecf477528b2f8e413881ab9a82aa2adc8471 Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Wed, 18 Jan 2023 12:55:42 +0100 Subject: [PATCH] Implement functions to compute and store mapping information for a vector of cells/faces --- include/deal.II/fe/mapping.h | 5 +- .../deal.II/matrix_free/fe_point_evaluation.h | 172 ++++- include/deal.II/non_matching/mapping_info.h | 677 +++++++++++++++++- tests/non_matching/mapping_info.cc | 477 ++++++++++++ tests/non_matching/mapping_info.output | 8 + 5 files changed, 1287 insertions(+), 52 deletions(-) create mode 100644 tests/non_matching/mapping_info.cc create mode 100644 tests/non_matching/mapping_info.output diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 9c2568abc6..a4b9578b15 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -55,7 +55,9 @@ namespace NonMatching { template class FEImmersedSurfaceValues; -} + template + class MappingInfo; +} // namespace NonMatching /** @@ -1305,6 +1307,7 @@ public: friend class FEFaceValues; friend class FESubfaceValues; friend class NonMatching::FEImmersedSurfaceValues; + friend class NonMatching::MappingInfo; }; diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 2d3aa93464..33232dba44 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -470,6 +470,20 @@ public: reinit(const typename Triangulation::cell_iterator &cell, const ArrayView> &unit_points); + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the cell in the MappingInfo object. + */ + void + reinit(const unsigned int cell_index); + + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the face in the MappingInfo object. + */ + void + reinit(const unsigned int cell_index, const unsigned int face_number); + /** * This function interpolates the finite element solution, represented by * `solution_values`, on the cell and `unit_points` passed to reinit(). @@ -585,6 +599,22 @@ public: DerivativeForm<1, spacedim, dim> inverse_jacobian(const unsigned int point_index) const; + /** + * Return the Jacobian determinant multiplied by the quadrature weight. This + * class or the MappingInfo object passed to this function needs to be + * constructed with UpdateFlags containing `update_JxW_values`. + */ + Number + JxW(const unsigned int point_index) const; + + /** + * Return the normal vector. This class or the MappingInfo object passed to + * this function needs to be constructed with UpdateFlags containing + * `update_normal_vectors`. + */ + Tensor<1, spacedim> + normal_vector(const unsigned int point_index) const; + /** * Return the position in real coordinates of the given point index among * the points passed to reinit(). @@ -599,6 +629,11 @@ public: Point unit_point(const unsigned int point_index) const; + /** + * Number of quadrature points of the current cell/face. + */ + const unsigned int n_q_points; + private: /** * Common setup function for both constructors. Does the setup for both fast @@ -713,10 +748,20 @@ private: */ SmartPointer> mapping_info; + /** + * The current cell index to access mapping data from mapping info. + */ + unsigned int current_cell_index; + + /** + * The current face number to access mapping data from mapping info. + */ + unsigned int current_face_number; + /** * The reference points specified at reinit(). */ - std::vector> unit_points; + ArrayView> unit_points; /** * Bool indicating if fast path is chosen. @@ -733,13 +778,16 @@ FEPointEvaluation::FEPointEvaluation( const FiniteElement &fe, const UpdateFlags update_flags, const unsigned int first_selected_component) - : mapping(&mapping) + : n_q_points(numbers::invalid_unsigned_int) + , mapping(&mapping) , fe(&fe) , update_flags(update_flags) , mapping_info_on_the_fly( std::make_unique>(mapping, update_flags)) , mapping_info(mapping_info_on_the_fly.get()) + , current_cell_index(numbers::invalid_unsigned_int) + , current_face_number(numbers::invalid_unsigned_int) { setup(first_selected_component); } @@ -751,10 +799,13 @@ FEPointEvaluation::FEPointEvaluation( NonMatching::MappingInfo &mapping_info, const FiniteElement & fe, const unsigned int first_selected_component) - : mapping(&mapping_info.get_mapping()) + : n_q_points(numbers::invalid_unsigned_int) + , mapping(&mapping_info.get_mapping()) , fe(&fe) , update_flags(mapping_info.get_update_flags()) , mapping_info(&mapping_info) + , current_cell_index(numbers::invalid_unsigned_int) + , current_face_number(numbers::invalid_unsigned_int) { setup(first_selected_component); } @@ -853,14 +904,45 @@ FEPointEvaluation::reinit( fe_values->reinit(cell); } - this->unit_points = - std::vector>(unit_points.begin(), unit_points.end()); + this->unit_points = unit_points; + + const_cast(n_q_points) = unit_points.size(); if (update_flags & update_values) - values.resize(unit_points.size(), numbers::signaling_nan()); + values.resize(n_q_points, numbers::signaling_nan()); if (update_flags & update_gradients) - gradients.resize(unit_points.size(), - numbers::signaling_nan()); + gradients.resize(n_q_points, numbers::signaling_nan()); +} + + + +template +void +FEPointEvaluation::reinit( + const unsigned int cell_index) +{ + current_cell_index = cell_index; + current_face_number = numbers::invalid_unsigned_int; + + const_cast(n_q_points) = + mapping_info->get_unit_points(current_cell_index, current_face_number) + .size(); +} + + + +template +void +FEPointEvaluation::reinit( + const unsigned int cell_index, + const unsigned int face_number) +{ + current_cell_index = cell_index; + current_face_number = face_number; + + const_cast(n_q_points) = + mapping_info->get_unit_points(current_cell_index, current_face_number) + .size(); } @@ -874,7 +956,8 @@ FEPointEvaluation::evaluate( const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr; if (precomputed_mapping) { - unit_points = mapping_info->get_unit_points(); + unit_points = + mapping_info->get_unit_points(current_cell_index, current_face_number); if (update_flags & update_values) values.resize(unit_points.size(), numbers::signaling_nan()); @@ -883,7 +966,7 @@ FEPointEvaluation::evaluate( numbers::signaling_nan()); } - if (unit_points.empty()) + if (unit_points.size() == 0) return; AssertDimension(solution_values.size(), fe->dofs_per_cell); @@ -946,11 +1029,12 @@ FEPointEvaluation::evaluate( Number>::set_gradient(val_and_grad.second, j, unit_gradients[i + j]); - gradients[i + j] = - apply_transformation(mapping_info->get_mapping_data() - .inverse_jacobians[i + j] - .transpose(), - unit_gradients[i + j]); + const auto &mapping_data = + mapping_info->get_mapping_data(current_cell_index, + current_face_number); + gradients[i + j] = apply_transformation( + mapping_data.inverse_jacobians[i + j].transpose(), + unit_gradients[i + j]); } } } @@ -1022,7 +1106,8 @@ FEPointEvaluation::integrate( const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr; if (precomputed_mapping) { - unit_points = mapping_info->get_unit_points(); + unit_points = + mapping_info->get_unit_points(current_cell_index, current_face_number); if (update_flags & update_values) values.resize(unit_points.size(), numbers::signaling_nan()); @@ -1086,9 +1171,12 @@ FEPointEvaluation::integrate( if (integration_flags & EvaluationFlags::gradients) for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) { - gradients[i + j] = apply_transformation( - mapping_info->get_mapping_data().inverse_jacobians[i + j], - gradients[i + j]); + const auto &mapping_data = + mapping_info->get_mapping_data(current_cell_index, + current_face_number); + gradients[i + j] = + apply_transformation(mapping_data.inverse_jacobians[i + j], + gradients[i + j]); internal::FEPointEvaluation:: EvaluatorTypeTraits::get_gradient( gradient, j, gradients[i + j]); @@ -1250,9 +1338,10 @@ inline DerivativeForm<1, dim, spacedim> FEPointEvaluation::jacobian( const unsigned int point_index) const { - AssertIndexRange(point_index, - mapping_info->get_mapping_data().jacobians.size()); - return mapping_info->get_mapping_data().jacobians[point_index]; + const auto &mapping_data = + mapping_info->get_mapping_data(current_cell_index, current_face_number); + AssertIndexRange(point_index, mapping_data.jacobians.size()); + return mapping_data.jacobians[point_index]; } @@ -1262,9 +1351,35 @@ inline DerivativeForm<1, spacedim, dim> FEPointEvaluation::inverse_jacobian( const unsigned int point_index) const { - AssertIndexRange(point_index, - mapping_info->get_mapping_data().inverse_jacobians.size()); - return mapping_info->get_mapping_data().inverse_jacobians[point_index]; + const auto &mapping_data = + mapping_info->get_mapping_data(current_cell_index, current_face_number); + AssertIndexRange(point_index, mapping_data.inverse_jacobians.size()); + return mapping_data.inverse_jacobians[point_index]; +} + + + +template +inline Number +FEPointEvaluation::JxW( + const unsigned int point_index) const +{ + const auto &mapping_data = + mapping_info->get_mapping_data(current_cell_index, current_face_number); + AssertIndexRange(point_index, mapping_data.JxW_values.size()); + return mapping_data.JxW_values[point_index]; +} + + +template +inline Tensor<1, spacedim> +FEPointEvaluation::normal_vector( + const unsigned int point_index) const +{ + const auto &mapping_data = + mapping_info->get_mapping_data(current_cell_index, current_face_number); + AssertIndexRange(point_index, mapping_data.normal_vectors.size()); + return mapping_data.normal_vectors[point_index]; } @@ -1274,9 +1389,10 @@ inline Point FEPointEvaluation::real_point( const unsigned int point_index) const { - AssertIndexRange(point_index, - mapping_info->get_mapping_data().quadrature_points.size()); - return mapping_info->get_mapping_data().quadrature_points[point_index]; + const auto &mapping_data = + mapping_info->get_mapping_data(current_cell_index, current_face_number); + AssertIndexRange(point_index, mapping_data.quadrature_points.size()); + return mapping_data.quadrature_points[point_index]; } diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index 241622d326..21347f0530 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -26,6 +26,7 @@ #include #include +#include #include #include #include @@ -45,6 +46,10 @@ namespace NonMatching class MappingInfo : public Subscriptor { public: + using MappingData = + dealii::internal::FEValuesImplementation::MappingRelatedData; + /** * Constructor. * @@ -57,27 +62,113 @@ namespace NonMatching MappingInfo(const Mapping &mapping, const UpdateFlags update_flags); /** - * Reinitialize the mapping information for the incoming cell and unit + * Compute the mapping information for the incoming cell and unit + * points. This overload is needed to resolve ambiguity. + */ + void + reinit(const typename Triangulation::cell_iterator &cell, + const std::vector> &unit_points); + + /** + * Compute the mapping information for the incoming cell and unit * points. */ void reinit(const typename Triangulation::cell_iterator &cell, const ArrayView> &unit_points); + /** + * Compute the mapping information for the given cell and + * quadrature formula. As opposed to the other `reinit` function, this + * method allows to access a `JxW` factor at the points. + */ + void + reinit(const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature); + + /** + * Compute the mapping information for the incoming vector of cells and + * corresponding vector of unit points. + * + * It is possible to give an IteratorRange to this + * function and together with @p n_unfiltered_cells specified this object + * can compress its storage while giving you access to the underlying data + * with the "unfiltered" index. The default argument for + * @p n_unfiltered_cells disables this built-in compression. + */ + template + void + reinit_cells( + const IteratorRange & cell_iterator_range, + const std::vector>> &unit_points_vector, + const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int); + + /** + * Compute the mapping information for the incoming vector of cells and + * corresponding vector of quadratures. As opposed to the other `reinit` + * function, this method allows to access a `JxW` factor at the points. + */ + template + void + reinit_cells( + const IteratorRange & cell_iterator_range, + const std::vector> &quadrature_vector, + const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int); + + /** + * Compute the mapping information for the incoming vector of cells and + * corresponding vector of ImmersedSurfaceQuadrature. + */ + template + void + reinit_surface( + const IteratorRange & cell_iterator_range, + const std::vector> &quadrature_vector, + const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int); + + /** + * Compute the mapping information for all faces of the incoming vector + * of cells and corresponding vector of quadratures. + */ + template + void + reinit_faces( + const IteratorRange & cell_iterator_range, + const std::vector>> &quadrature_vector, + const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int); + /** * Getter function for current unit points. + * + * @p cell_index and @p face_number are the indices + * into the compressed, CRS like data storage of unit points. + * + * If you have initialized this object with reinit_cells() you can access + * the stored unit points of the cell with the respective @p cell_index + * (and the default argument for @p face_number). + * + * If you have initialized this object with reinit_faces() you can access + * the stored unit points of the faces on the cell with the respective @p cell_index + * and the respective local @p face_number. + * + * If you have initialized this object with reinit() you can access the + * stored unit points of a single cell with the default arguments. + * + * The correct state of this object is checked in this call (in debug mode). */ - const std::vector> & - get_unit_points() const; + const ArrayView> + get_unit_points( + const unsigned int cell_index = numbers::invalid_unsigned_int, + const unsigned int face_number = numbers::invalid_unsigned_int) const; /** * Getter function for computed mapping data. This function accesses * internal data and is therefore not a stable interface. */ - const dealii::internal::FEValuesImplementation::MappingRelatedData - & - get_mapping_data() const; + const MappingData & + get_mapping_data( + const unsigned int cell_index = numbers::invalid_unsigned_int, + const unsigned int face_number = numbers::invalid_unsigned_int) const; /** * Getter function for underlying mapping. @@ -92,6 +183,23 @@ namespace NonMatching get_update_flags() const; private: + /** + * Enum class for reinitialized states. + */ + enum class State + { + invalid, + single_cell, + cell_vector, + faces_on_cells_in_vector + }; + + /** + * Enum class that stores the currently initialized state + * upon the last call of reinit(). + */ + State state; + /** * Compute the mapping related data for the given @p mapping, * @p cell and @p unit_points that is required by the FEPointEvaluation @@ -100,13 +208,42 @@ namespace NonMatching void compute_mapping_data_for_generic_points( const typename Triangulation::cell_iterator &cell, - const ArrayView> & unit_points); + const ArrayView> & unit_points, + MappingData & mapping_data); + + /** + * Compute the mapping related data for the given @p mapping, + * @p cell and @p quadrature that is required by the FEPointEvaluation + * class. + */ + void + compute_mapping_data_for_immersed_surface_quadrature( + const typename Triangulation::cell_iterator &cell, + const ImmersedSurfaceQuadrature & quadrature, + MappingData & mapping_data); + + /** + * Compute the mapping related data for the given @p mapping, @p cell, + * @p face_no and @p quadrature that is required by the FEPointEvaluation + * class. + */ + void + compute_mapping_data_for_face_quadrature( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & quadrature, + MappingData & mapping_data); /** * The reference points specified at reinit(). */ std::vector> unit_points; + /** + * Offset to point to the first unit point of a cell/face + */ + std::vector unit_points_index; + /** * A pointer to the underlying mapping. */ @@ -126,8 +263,24 @@ namespace NonMatching * The internal data container for mapping information. The implementation * is subject to future changes. */ - dealii::internal::FEValuesImplementation::MappingRelatedData - mapping_data; + std::vector mapping_data; + + /** + * Offset to point to the first element of a cell in internal data + * containers. + */ + std::vector cell_index_offset; + + /** + * A vector that converts the cell index to a compressed cell index for e.g. + * a filtered IteratorRange. + */ + std::vector cell_index_to_compressed_cell_index; + + /** + * A bool that determines weather cell index compression should be done. + */ + bool do_cell_index_compression; }; // ----------------------- template functions ---------------------- @@ -141,8 +294,12 @@ namespace NonMatching { update_flags_mapping = update_default; // translate update flags - if (update_flags & update_jacobians) + if (update_flags & update_jacobians || update_flags & update_JxW_values) update_flags_mapping |= update_jacobians; + if (update_flags & update_JxW_values) + update_flags_mapping |= update_JxW_values; + if (update_flags & update_normal_vectors) + update_flags_mapping |= update_normal_vectors; if (update_flags & update_gradients || update_flags & update_inverse_jacobians) update_flags_mapping |= update_inverse_jacobians; @@ -157,30 +314,451 @@ namespace NonMatching void MappingInfo::reinit( const typename Triangulation::cell_iterator &cell, - const ArrayView> & unit_points) + const std::vector> & unit_points_in) + { + reinit(cell, make_array_view(unit_points_in.begin(), unit_points_in.end())); + } + + + + template + void + MappingInfo::reinit( + 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; + } + + + + template + void + MappingInfo::reinit( + const typename Triangulation::cell_iterator &cell, + const Quadrature & quadrature) + { + const auto &points = quadrature.get_points(); + const auto &weights = quadrature.get_weights(); + + reinit(cell, points); + + 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]; + } + + + + template + template + void + MappingInfo::reinit_cells( + const IteratorRange & cell_iterator_range, + const std::vector>> &unit_points_vector, + const unsigned int n_unfiltered_cells) { - this->unit_points = - std::vector>(unit_points.begin(), unit_points.end()); - compute_mapping_data_for_generic_points(cell, unit_points); + do_cell_index_compression = + n_unfiltered_cells != numbers::invalid_unsigned_int; + + const unsigned int n_cells = unit_points_vector.size(); + AssertDimension(n_cells, + std::distance(cell_iterator_range.begin(), + cell_iterator_range.end())); + + // 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()); + + const unsigned int n_unit_points = unit_points_index.back(); + + unit_points.resize(n_unit_points); + mapping_data.resize(n_cells); + + if (do_cell_index_compression) + cell_index_to_compressed_cell_index.resize(n_unfiltered_cells, + numbers::invalid_unsigned_int); + 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]) + { + *it = unit_point; + ++it; + } + + compute_mapping_data_for_generic_points(cell, + unit_points_vector[cell_index], + mapping_data[cell_index]); + + if (do_cell_index_compression) + cell_index_to_compressed_cell_index[cell->active_cell_index()] = + cell_index; + + ++cell_index; + } + + state = State::cell_vector; } template - const std::vector> & - MappingInfo::get_unit_points() const + template + void + MappingInfo::reinit_cells( + const IteratorRange & cell_iterator_range, + const std::vector> &quadrature_vector, + const unsigned int n_unfiltered_cells) { - return unit_points; + 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 - const dealii::internal::FEValuesImplementation::MappingRelatedData & - MappingInfo::get_mapping_data() const + template + void + MappingInfo::reinit_surface( + const IteratorRange & cell_iterator_range, + const std::vector> &quadrature_vector, + const unsigned int n_unfiltered_cells) { - return mapping_data; + do_cell_index_compression = + n_unfiltered_cells != numbers::invalid_unsigned_int; + + if (update_flags_mapping & (update_JxW_values | update_normal_vectors)) + update_flags_mapping |= update_covariant_transformation; + + const unsigned int n_cells = quadrature_vector.size(); + AssertDimension(n_cells, + std::distance(cell_iterator_range.begin(), + cell_iterator_range.end())); + + // fill unit points index offset vector + unit_points_index.reserve(n_cells + 1); + unit_points_index.push_back(0); + for (const auto &quadrature : quadrature_vector) + unit_points_index.push_back(unit_points_index.back() + + quadrature.get_points().size()); + + const unsigned int n_unit_points = unit_points_index.back(); + + unit_points.resize(n_unit_points); + mapping_data.resize(n_cells); + + if (do_cell_index_compression) + cell_index_to_compressed_cell_index.resize(n_unfiltered_cells, + numbers::invalid_unsigned_int); + unsigned int cell_index = 0; + for (const auto &cell : cell_iterator_range) + { + const auto &quadrature = quadrature_vector[cell_index]; + + auto it = unit_points.begin() + unit_points_index[cell_index]; + for (const auto &unit_point : quadrature.get_points()) + { + *it = unit_point; + ++it; + } + + compute_mapping_data_for_immersed_surface_quadrature( + cell, quadrature, mapping_data[cell_index]); + + if (do_cell_index_compression) + cell_index_to_compressed_cell_index[cell->active_cell_index()] = + cell_index; + + ++cell_index; + } + + state = State::cell_vector; + } + + + + template + template + void + MappingInfo::reinit_faces( + 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 = quadrature_vector.size(); + AssertDimension(n_cells, + std::distance(cell_iterator_range.begin(), + cell_iterator_range.end())); + + // fill cell index offset vector + cell_index_offset.resize(n_cells); + unsigned int n_faces = 0; + unsigned int cell_index = 0; + for (const auto &cell : cell_iterator_range) + { + cell_index_offset[cell_index] = n_faces; + n_faces += cell->n_faces(); + ++cell_index; + } + + // fill unit points index offset vector + unit_points_index.resize(n_faces + 1); + cell_index = 0; + unsigned int n_unit_points = 0; + for (const auto &cell : cell_iterator_range) + { + for (const auto &f : cell->face_indices()) + { + const unsigned int current_face_index = + cell_index_offset[cell_index] + f; + + unit_points_index[current_face_index] = n_unit_points; + n_unit_points += + quadrature_vector[cell_index][f].get_points().size(); + } + + ++cell_index; + } + unit_points_index[n_faces] = n_unit_points; + + // compress indices + if (do_cell_index_compression) + cell_index_to_compressed_cell_index.resize(n_unfiltered_cells, + numbers::invalid_unsigned_int); + + // fill unit points and mapping data for every face of all cells + unit_points.resize(n_unit_points); + mapping_data.resize(n_faces); + cell_index = 0; + QProjector q_projector; + for (const auto &cell : cell_iterator_range) + { + const auto &quadratures_on_faces = quadrature_vector[cell_index]; + + Assert(quadratures_on_faces.size() == cell->n_faces(), + ExcDimensionMismatch(quadratures_on_faces.size(), + cell->n_faces())); + + for (const auto &f : cell->face_indices()) + { + const auto &quadrature_on_face = quadratures_on_faces[f]; + + const auto quadrature_on_cell = + q_projector.project_to_face(cell->reference_cell(), + quadrature_on_face, + f); + + const auto &unit_points_on_cell = quadrature_on_cell.get_points(); + + const unsigned int current_face_index = + cell_index_offset[cell_index] + f; + + auto it = + unit_points.begin() + unit_points_index[current_face_index]; + for (const auto &unit_point : unit_points_on_cell) + { + *it = unit_point; + ++it; + } + + compute_mapping_data_for_face_quadrature( + cell, f, quadrature_on_face, mapping_data[current_face_index]); + } + if (do_cell_index_compression) + cell_index_to_compressed_cell_index[cell->active_cell_index()] = + cell_index; + + ++cell_index; + } + + state = State::faces_on_cells_in_vector; + } + + + + template + const ArrayView> + MappingInfo::get_unit_points( + const unsigned int cell_index, + const unsigned int face_number) const + { + if (cell_index == numbers::invalid_unsigned_int && + face_number == numbers::invalid_unsigned_int) + { + Assert(state == State::single_cell, + ExcMessage( + "This mapping info is not reinitialized for a single cell!")); + return unit_points; + } + else if (face_number == numbers::invalid_unsigned_int) + { + Assert(state == State::cell_vector, + ExcMessage( + "This mapping info is not reinitialized for a cell vector!")); + if (do_cell_index_compression) + { + Assert(cell_index_to_compressed_cell_index[cell_index] != + numbers::invalid_unsigned_int, + ExcMessage("Mapping info object was not initialized for this" + " active cell index!")); + const auto it_begin = + unit_points.begin() + + unit_points_index + [cell_index_to_compressed_cell_index[cell_index]]; + const auto it_end = + unit_points.begin() + + unit_points_index + [cell_index_to_compressed_cell_index[cell_index] + 1]; + return make_array_view(it_begin, it_end); + } + else + { + const auto it_begin = + unit_points.begin() + unit_points_index[cell_index]; + const auto it_end = + unit_points.begin() + unit_points_index[cell_index + 1]; + return make_array_view(it_begin, it_end); + } + } + else if (cell_index != numbers::invalid_unsigned_int) + { + Assert(state == State::faces_on_cells_in_vector, + ExcMessage("This mapping info is not reinitialized for faces" + " on cells in a vector!")); + if (do_cell_index_compression) + { + Assert( + cell_index_to_compressed_cell_index[cell_index] != + numbers::invalid_unsigned_int, + ExcMessage( + "Mapping info object was not initialized for this active cell index" + " and corresponding face numbers!")); + const unsigned int current_face_index = + cell_index_offset + [cell_index_to_compressed_cell_index[cell_index]] + + face_number; + const auto it_begin = + unit_points.begin() + unit_points_index[current_face_index]; + const auto it_end = + unit_points.begin() + unit_points_index[current_face_index + 1]; + return make_array_view(it_begin, it_end); + } + else + { + const unsigned int current_face_index = + cell_index_offset[cell_index] + face_number; + const auto it_begin = + unit_points.begin() + unit_points_index[current_face_index]; + const auto it_end = + unit_points.begin() + unit_points_index[current_face_index + 1]; + return make_array_view(it_begin, it_end); + } + } + else + AssertThrow( + false, + ExcMessage( + "cell_index has to be specified if face_number is specified!")); + } + + + + template + const typename MappingInfo::MappingData & + MappingInfo::get_mapping_data( + const unsigned int cell_index, + const unsigned int face_number) const + { + if (cell_index == numbers::invalid_unsigned_int && + face_number == numbers::invalid_unsigned_int) + { + Assert(state == State::single_cell, + ExcMessage( + "This mapping info is not reinitialized for a single cell!")); + return mapping_data[0]; + } + else if (face_number == numbers::invalid_unsigned_int) + { + Assert(state == State::cell_vector, + ExcMessage( + "This mapping info is not reinitialized for a cell vector!")); + if (do_cell_index_compression) + { + Assert(cell_index_to_compressed_cell_index[cell_index] != + numbers::invalid_unsigned_int, + ExcMessage("Mapping info object was not initialized for this" + " active cell index!")); + return mapping_data + [cell_index_to_compressed_cell_index[cell_index]]; + } + else + return mapping_data[cell_index]; + } + else if (cell_index != numbers::invalid_unsigned_int) + { + Assert(state == State::faces_on_cells_in_vector, + ExcMessage("This mapping info is not reinitialized for faces" + " on cells in a vector!")); + if (do_cell_index_compression) + { + Assert( + cell_index_to_compressed_cell_index[cell_index] != + numbers::invalid_unsigned_int, + ExcMessage( + "Mapping info object was not initialized for this active cell index" + " and corresponding face numbers!")); + return mapping_data + [cell_index_offset + [cell_index_to_compressed_cell_index[cell_index]] + + face_number]; + } + else + return mapping_data[cell_index_offset[cell_index] + face_number]; + } + else + AssertThrow( + false, + ExcMessage( + "cell_index has to be specified if face_number is specified!")); } @@ -207,7 +785,8 @@ namespace NonMatching void MappingInfo::compute_mapping_data_for_generic_points( const typename Triangulation::cell_iterator &cell, - const ArrayView> & unit_points) + const ArrayView> & unit_points, + MappingData & mapping_data) { if (const MappingQ *mapping_q = dynamic_cast *>(&(*mapping))) @@ -246,6 +825,58 @@ namespace NonMatching mapping_data.quadrature_points[q] = fe_values.quadrature_point(q); } } + + + + template + void + MappingInfo:: + compute_mapping_data_for_immersed_surface_quadrature( + const typename Triangulation::cell_iterator &cell, + const ImmersedSurfaceQuadrature & quadrature, + MappingData & mapping_data) + { + 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_data(update_flags_mapping, quadrature); + + mapping->fill_fe_immersed_surface_values(cell, + quadrature, + *internal_mapping_data, + mapping_data); + } + + + + template + void + MappingInfo::compute_mapping_data_for_face_quadrature( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & quadrature, + MappingData & mapping_data) + { + 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->fill_fe_face_values(cell, + face_no, + hp::QCollection(quadrature), + *internal_mapping_data, + mapping_data); + } } // namespace NonMatching DEAL_II_NAMESPACE_CLOSE diff --git a/tests/non_matching/mapping_info.cc b/tests/non_matching/mapping_info.cc new file mode 100644 index 0000000000..f6c7988fea --- /dev/null +++ b/tests/non_matching/mapping_info.cc @@ -0,0 +1,477 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the NonMatching::MappingInfo class together with FEPointEvaluation and + * compare to NonMatching::FEValues + */ + +#include + +#include + +#include + +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +void +test(const bool filtered_compression) +{ + constexpr unsigned int dim = 2; + constexpr unsigned int degree = 1; + + FE_Q fe_q(degree); + + Triangulation tria; + + MappingQ mapping(degree); + + GridGenerator::subdivided_hyper_cube(tria, 4); + + DoFHandler dof_handler(tria); + + dof_handler.distribute_dofs(fe_q); + + Functions::SignedDistance::Sphere level_set; + + Vector level_set_vec(dof_handler.n_dofs()); + + VectorTools::interpolate(dof_handler, level_set, level_set_vec); + + NonMatching::MeshClassifier mesh_classifier(dof_handler, level_set_vec); + mesh_classifier.reclassify(); + + hp::QCollection<1> q_collection((QGauss<1>(degree))); + + NonMatching::DiscreteQuadratureGenerator quadrature_generator( + q_collection, dof_handler, level_set_vec); + + NonMatching::DiscreteFaceQuadratureGenerator face_quadrature_generator( + q_collection, dof_handler, level_set_vec); + + // FEPointEvaluation + NonMatching::MappingInfo mapping_info_cell( + mapping, update_values | update_gradients | update_JxW_values); + + NonMatching::MappingInfo mapping_info_surface(mapping, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors); + + NonMatching::MappingInfo mapping_info_faces(mapping, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors); + + auto physical = + [&](const typename Triangulation::active_cell_iterator &i) { + return (mesh_classifier.location_to_level_set(i) == + NonMatching::LocationToLevelSet::intersected) || + (mesh_classifier.location_to_level_set(i) == + NonMatching::LocationToLevelSet::inside); + }; + + if (filtered_compression) + { + const auto physical_active_cell_iterators = + tria.active_cell_iterators() | physical; + + unsigned int n_physical_cells = 0; + for (const auto &cell : physical_active_cell_iterators) + { + ++n_physical_cells; + } + + std::vector> quad_vec_cell; + quad_vec_cell.reserve(n_physical_cells); + std::vector> quad_vec_surface; + quad_vec_surface.reserve(n_physical_cells); + std::vector>> quad_vec_faces( + n_physical_cells); + + unsigned int counter = 0; + for (const auto &cell : physical_active_cell_iterators) + { + quadrature_generator.generate(cell); + quad_vec_cell.push_back(quadrature_generator.get_inside_quadrature()); + quad_vec_surface.push_back( + quadrature_generator.get_surface_quadrature()); + + for (auto f : cell->face_indices()) + { + face_quadrature_generator.generate(cell, f); + quad_vec_faces[counter].push_back( + face_quadrature_generator.get_inside_quadrature()); + } + + ++counter; + } + + mapping_info_cell.reinit_cells(physical_active_cell_iterators, + quad_vec_cell, + tria.n_active_cells()); + mapping_info_surface.reinit_surface(physical_active_cell_iterators, + quad_vec_surface, + tria.n_active_cells()); + mapping_info_faces.reinit_faces(physical_active_cell_iterators, + quad_vec_faces, + tria.n_active_cells()); + } + else + { + std::vector> quad_vec_cell; + std::vector> quad_vec_surface; + std::vector>> quad_vec_faces( + tria.n_active_cells()); + for (const auto &cell : tria.active_cell_iterators()) + { + quadrature_generator.generate(cell); + quad_vec_cell.push_back(quadrature_generator.get_inside_quadrature()); + quad_vec_surface.push_back( + quadrature_generator.get_surface_quadrature()); + + for (auto f : cell->face_indices()) + { + face_quadrature_generator.generate(cell, f); + quad_vec_faces[cell->active_cell_index()].push_back( + face_quadrature_generator.get_inside_quadrature()); + } + } + + mapping_info_cell.reinit_cells(tria.active_cell_iterators(), + quad_vec_cell); + mapping_info_surface.reinit_surface(tria.active_cell_iterators(), + quad_vec_surface); + mapping_info_faces.reinit_faces(tria.active_cell_iterators(), + quad_vec_faces); + } + + Vector src(dof_handler.n_dofs()), dst_cell(dof_handler.n_dofs()), + dst_surface(dof_handler.n_dofs()), dst_faces(dof_handler.n_dofs()); + + for (auto &v : src) + v = random_value(); + + FEPointEvaluation<1, dim, dim, double> fe_point_cell(mapping_info_cell, fe_q); + FEPointEvaluation<1, dim, dim, double> fe_point_surface(mapping_info_surface, + fe_q); + FEPointEvaluation<1, dim, dim, double> fe_point_faces_m(mapping_info_faces, + fe_q); + FEPointEvaluation<1, dim, dim, double> fe_point_faces_p(mapping_info_faces, + fe_q); + + std::vector solution_values_in(fe_q.dofs_per_cell); + std::vector solution_values_neighbor_in(fe_q.dofs_per_cell); + std::vector solution_values_cell_out(fe_q.dofs_per_cell); + std::vector solution_values_surface_out(fe_q.dofs_per_cell); + std::vector solution_values_faces_out(fe_q.dofs_per_cell); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (!(physical(cell))) + continue; + + cell->get_dof_values(src, + solution_values_in.begin(), + solution_values_in.end()); + + auto test_fe_point = [&](FEPointEvaluation<1, dim, dim, double> &fe_point, + std::vector &solution_values_out) { + fe_point.reinit(cell->active_cell_index()); + + fe_point.evaluate(solution_values_in, + EvaluationFlags::values | EvaluationFlags::gradients); + + for (unsigned int q = 0; q < fe_point.n_q_points; ++q) + { + fe_point.submit_value(fe_point.JxW(q) * fe_point.get_value(q), q); + fe_point.submit_gradient(fe_point.JxW(q) * fe_point.get_gradient(q), + q); + } + + fe_point.integrate(solution_values_out, + EvaluationFlags::values | + EvaluationFlags::gradients); + }; + + test_fe_point(fe_point_cell, solution_values_cell_out); + test_fe_point(fe_point_surface, solution_values_surface_out); + + for (const auto f : cell->face_indices()) + { + if (cell->at_boundary(f) || !physical(cell->neighbor(f))) + continue; + + fe_point_faces_m.reinit(cell->active_cell_index(), f); + fe_point_faces_p.reinit(cell->neighbor(f)->active_cell_index(), + cell->neighbor_of_neighbor(f)); + + cell->neighbor(f)->get_dof_values(src, + solution_values_neighbor_in.begin(), + solution_values_neighbor_in.end()); + + fe_point_faces_m.evaluate(solution_values_in, + EvaluationFlags::values | + EvaluationFlags::gradients); + + fe_point_faces_p.evaluate(solution_values_neighbor_in, + EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < fe_point_faces_m.n_q_points; ++q) + { + fe_point_faces_m.submit_value(fe_point_faces_m.JxW(q) * + (fe_point_faces_m.get_value(q) - + fe_point_faces_p.get_value(q)), + q); + fe_point_faces_m.submit_gradient( + fe_point_faces_m.JxW(q) * (fe_point_faces_m.get_gradient(q) - + fe_point_faces_p.get_gradient(q)), + q); + } + + fe_point_faces_m.integrate(solution_values_faces_out, + EvaluationFlags::values | + EvaluationFlags::gradients); + + cell->distribute_local_to_global( + Vector(solution_values_faces_out.begin(), + solution_values_faces_out.end()), + dst_faces); + } + + cell->distribute_local_to_global( + Vector(solution_values_cell_out.begin(), + solution_values_cell_out.end()), + dst_cell); + + cell->distribute_local_to_global( + Vector(solution_values_surface_out.begin(), + solution_values_surface_out.end()), + dst_surface); + } + + + // FEValues + const QGauss<1> quadrature_1D(degree); + + NonMatching::RegionUpdateFlags region_update_flags; + region_update_flags.inside = update_values | update_gradients | + update_JxW_values | update_quadrature_points; + region_update_flags.surface = update_values | update_gradients | + update_JxW_values | update_quadrature_points | + update_normal_vectors; + + hp::FECollection fe_collection(fe_q); + + NonMatching::FEValues non_matching_fe_values(fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + dof_handler, + level_set_vec); + + NonMatching::FEInterfaceValues non_matching_fe_interface_values( + fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + dof_handler, + level_set_vec); + + Vector dst_cell_2(dof_handler.n_dofs()), + dst_surface_2(dof_handler.n_dofs()), dst_faces_2(dof_handler.n_dofs()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + non_matching_fe_values.reinit(cell); + + cell->get_dof_values(src, + solution_values_in.begin(), + solution_values_in.end()); + + const auto &inside_fe_values = + non_matching_fe_values.get_inside_fe_values(); + + const auto &surface_fe_values = + non_matching_fe_values.get_surface_fe_values(); + + + auto test_fe_values = [&](const auto & fe_values, + std::vector &solution_values_out, + Vector & dst) { + if (fe_values) + { + std::vector> solution_gradients( + fe_values->n_quadrature_points); + std::vector solution_values(fe_values->n_quadrature_points); + + for (const auto q : fe_values->quadrature_point_indices()) + { + double values = 0.; + Tensor<1, dim> gradients; + + for (const auto i : fe_values->dof_indices()) + { + gradients += + solution_values_in[i] * fe_values->shape_grad(i, q); + values += + solution_values_in[i] * fe_values->shape_value(i, q); + } + solution_gradients[q] = gradients * fe_values->JxW(q); + solution_values[q] = values * fe_values->JxW(q); + } + + for (const auto i : fe_values->dof_indices()) + { + double sum_gradients = 0.; + double sum_values = 0.; + for (const auto q : fe_values->quadrature_point_indices()) + { + sum_gradients += + solution_gradients[q] * fe_values->shape_grad(i, q); + sum_values += + solution_values[q] * fe_values->shape_value(i, q); + } + + solution_values_cell_out[i] = sum_gradients + sum_values; + } + + cell->distribute_local_to_global( + Vector(solution_values_cell_out.begin(), + solution_values_cell_out.end()), + dst); + } + }; + + test_fe_values(inside_fe_values, solution_values_cell_out, dst_cell_2); + test_fe_values(surface_fe_values, + solution_values_surface_out, + dst_surface_2); + + for (const auto f : cell->face_indices()) + { + if (cell->at_boundary(f)) + continue; + + non_matching_fe_interface_values.reinit( + cell, + f, + numbers::invalid_unsigned_int, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + numbers::invalid_unsigned_int); + + const auto &fe_interface_values = + non_matching_fe_interface_values.get_inside_fe_values(); + + if (fe_interface_values) + { + const auto &fe_face_values_m = + fe_interface_values->get_fe_face_values(0); + + const auto &fe_face_values_p = + fe_interface_values->get_fe_face_values(1); + + cell->neighbor(f)->get_dof_values( + src, + solution_values_neighbor_in.begin(), + solution_values_neighbor_in.end()); + + std::vector> solution_gradients( + fe_face_values_m.n_quadrature_points); + std::vector solution_values( + fe_face_values_m.n_quadrature_points); + + for (const auto q : fe_face_values_m.quadrature_point_indices()) + { + double values = 0.; + Tensor<1, dim> gradients; + + for (const auto i : fe_face_values_m.dof_indices()) + { + gradients += solution_values_in[i] * + fe_face_values_m.shape_grad(i, q) - + solution_values_neighbor_in[i] * + fe_face_values_p.shape_grad(i, q); + values += solution_values_in[i] * + fe_face_values_m.shape_value(i, q) - + solution_values_neighbor_in[i] * + fe_face_values_p.shape_value(i, q); + } + solution_gradients[q] = gradients * fe_face_values_m.JxW(q); + solution_values[q] = values * fe_face_values_m.JxW(q); + } + + for (const auto i : fe_face_values_m.dof_indices()) + { + double sum_gradients = 0.; + double sum_values = 0.; + for (const auto q : + fe_face_values_m.quadrature_point_indices()) + { + sum_gradients += solution_gradients[q] * + fe_face_values_m.shape_grad(i, q); + sum_values += + solution_values[q] * fe_face_values_m.shape_value(i, q); + } + + solution_values_faces_out[i] = sum_gradients + sum_values; + } + + cell->distribute_local_to_global( + Vector(solution_values_faces_out.begin(), + solution_values_faces_out.end()), + dst_faces_2); + } + } + } + + deallog << "check difference l2 norm cell: " + << dst_cell.l2_norm() - dst_cell_2.l2_norm() << std::endl; + + deallog << "check difference l2 norm surface: " + << dst_surface.l2_norm() - dst_surface_2.l2_norm() << std::endl; + + deallog << "check difference l2 norm faces: " + << dst_faces.l2_norm() - dst_faces_2.l2_norm() << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + + initlog(); + + test(true); + deallog << std::endl; + test(false); +} diff --git a/tests/non_matching/mapping_info.output b/tests/non_matching/mapping_info.output new file mode 100644 index 0000000000..b4737523f9 --- /dev/null +++ b/tests/non_matching/mapping_info.output @@ -0,0 +1,8 @@ + +DEAL::check difference l2 norm cell: 0.00000 +DEAL::check difference l2 norm surface: 0.00000 +DEAL::check difference l2 norm faces: 0.00000 +DEAL:: +DEAL::check difference l2 norm cell: 0.00000 +DEAL::check difference l2 norm surface: 0.00000 +DEAL::check difference l2 norm faces: 0.00000 -- 2.39.5