From: Peter Munch Date: Sat, 5 Aug 2023 12:44:42 +0000 (+0200) Subject: FEPointEvaluation/MappingInfo for simplices X-Git-Tag: relicensing~559^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=61c1e0f93b52bbee0686db884aaa7796c3ffdee9;p=dealii.git FEPointEvaluation/MappingInfo for simplices --- diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 4f71313c46..841adb4725 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -1603,6 +1603,32 @@ FEPointEvaluation::reinit( current_face_number = numbers::invalid_unsigned_int; do_reinit(); + + if (!fast_path) + { + const auto unit_points_vectorized = mapping_info->get_unit_point( + mapping_info->compute_unit_point_index_offset(current_cell_index, + current_face_number)); + const unsigned int n_q_points_unvectorized = + mapping_info->get_n_q_points_unvectorized(current_cell_index, + current_face_number); + + std::vector> unit_points(n_q_points_unvectorized); + + for (unsigned int v = 0; v < n_q_points_unvectorized; ++v) + for (unsigned int d = 0; d < dim; ++d) + unit_points[v][d] = unit_points_vectorized[v / n_lanes_internal][d] + [v % n_lanes_internal]; + + fe_values = std::make_shared>( + *mapping, + *fe, + Quadrature( + std::vector>(unit_points.begin(), unit_points.end())), + update_flags); + + fe_values->reinit(mapping_info->get_cell_iterator(current_cell_index)); + } } diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index 532dbe635b..b8c571fd95 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -216,8 +216,10 @@ namespace NonMatching /** * Constructor which sets the default arguments. */ - AdditionalData(const bool use_global_weights = false) + AdditionalData(const bool use_global_weights = false, + const bool store_cells = false) : use_global_weights(use_global_weights) + , store_cells(store_cells) {} /** @@ -226,6 +228,15 @@ namespace NonMatching * QSimplex::compute_affine_transformation(). */ bool use_global_weights; + + /** + * During the reinit() function calls, cells are passed as + * argument. In the default case, the cell is not stored, + * since all relevant mapping related information is precomputed. + * Hoever, this flag enables that the cells are stored so that + * they can be accessed later on. + */ + bool store_cells; }; /** @@ -422,6 +433,14 @@ namespace NonMatching get_n_q_points_unvectorized(const unsigned int cell_index, const unsigned int face_number) const; + /** + * Return cell iterator. + * + * @note This call is only possible if AdditionalData::store_cells is enabled. + */ + typename Triangulation::cell_iterator + get_cell_iterator(const unsigned int cell_index) const; + private: using MappingData = dealii::internal::FEValuesImplementation::MappingRelatedData is_reinitialized; + + /** + * Reference to the triangulation passed via the cells to the + * reinit functions. This field is only set if + * AdditionalData::store_cells is enabled. + */ + SmartPointer> triangulation; + + /** + * Level and indices of cells passed to the reinit functions. This + * vector is only filled if AdditionalData::store_cells is enabled. + */ + std::vector> cell_level_and_indices; }; // ----------------------- template functions ---------------------- @@ -807,6 +839,9 @@ namespace NonMatching n_q_points_unvectorized.reserve(n_cells); + if (additional_data.store_cells) + cell_level_and_indices.resize(n_cells); + // fill unit points index offset vector unit_points_index.reserve(n_cells + 1); unit_points_index.push_back(0); @@ -844,6 +879,12 @@ namespace NonMatching unsigned int cell_index = 0; for (const auto &cell : cell_iterator_range) { + if (additional_data.store_cells) + { + this->triangulation = &cell->get_triangulation(); + cell_level_and_indices[cell_index] = {cell->level(), cell->index()}; + } + // store unit points const unsigned int n_q_points = compute_n_q_points( n_q_points_unvectorized[cell_index]); @@ -1164,6 +1205,22 @@ namespace NonMatching + template + typename Triangulation::cell_iterator + MappingInfo::get_cell_iterator( + const unsigned int cell_index) const + { + Assert( + additional_data.store_cells, + ExcMessage( + "Cells have been not stored. You can enable this by Additional::store_cells.")); + return {triangulation.get(), + cell_level_and_indices[cell_index].first, + cell_level_and_indices[cell_index].second}; + } + + + template template unsigned int