From af5aa8a2fdbaf24624e62615de73a46ae4d98353 Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Thu, 18 Jan 2024 22:18:29 +0100 Subject: [PATCH] Address review comments --- ...aluation_template_face_factory.templates.h | 1 + include/deal.II/matrix_free/fe_evaluation.h | 65 +- .../deal.II/matrix_free/fe_point_evaluation.h | 2291 +++++++++-------- include/deal.II/non_matching/mapping_info.h | 6 +- tests/non_matching/mapping_info_04.cc | 4 +- tests/non_matching/mapping_info_05.cc | 6 +- 6 files changed, 1216 insertions(+), 1157 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h index 79519e39a8..608d03e1b6 100644 --- a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h @@ -47,6 +47,7 @@ namespace internal } + template void FEFaceEvaluationFactory::project_to_face( diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 1ed360ce38..dfe94b4dd4 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2750,20 +2750,18 @@ public: boundary_id() const; /** - * The number of degrees of freedom of a single component on the cell for - * the underlying evaluation object. Usually close to - * static_dofs_per_component, but the number depends on the actual element - * selected and is thus not static. + * Get the number of degrees of freedom of a single component which are + * projected onto a face. */ - const unsigned int dofs_per_component_on_cell; + unsigned int + get_dofs_per_component_projected_to_face(); /** - * The number of degrees of freedom on the cell accumulated over all - * components in the current evaluation object. Usually close to - * static_dofs_per_cell = static_dofs_per_component*n_components, but the - * number depends on the actual element selected and is thus not static. + * Get the number of degrees of freedom accumulated over all + * components which are projected onto a face. */ - const unsigned int dofs_per_cell; + unsigned int + get_dofs_projected_to_face(); /** * The number of degrees of freedom of a single component on the cell for @@ -2771,7 +2769,7 @@ public: * static_dofs_per_component, but the number depends on the actual element * selected and is thus not static. */ - const unsigned int dofs_per_component_on_face; + const unsigned int dofs_per_component; /** * The number of degrees of freedom on the cell accumulated over all @@ -2779,7 +2777,7 @@ public: * static_dofs_per_cell = static_dofs_per_component*n_components, but the * number depends on the actual element selected and is thus not static. */ - const unsigned int dofs_per_face; + const unsigned int dofs_per_cell; /** * The number of quadrature points in use. If the number of quadrature @@ -8310,10 +8308,8 @@ inline FEFaceEvaluationdata->dofs_per_component_on_cell) + , dofs_per_component(this->data->dofs_per_component_on_cell) , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_) - , dofs_per_component_on_face(this->data->dofs_per_component_on_face) - , dofs_per_face(this->data->dofs_per_component_on_face * n_components_) , n_q_points(this->n_quadrature_points) {} @@ -9414,6 +9410,45 @@ FEFaceEvaluation +unsigned int +FEFaceEvaluation< + dim, + fe_degree, + n_q_points_1d, + n_components_, + Number, + VectorizedArrayType>::get_dofs_per_component_projected_to_face() +{ + return this->data->dofs_per_component_on_face; +} + + + +template +unsigned int +FEFaceEvaluation::get_dofs_projected_to_face() +{ + return this->data->dofs_per_component_on_face * n_components_; +} + + + /*------------------------- end FEFaceEvaluation ------------------------- */ diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index da845ef1bd..8ff84419aa 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -1111,1004 +1111,1309 @@ protected: const bool is_interior; }; -// ----------------------- template and inline function ---------------------- -template -FEPointEvaluationBase:: - FEPointEvaluationBase(const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component) - : n_q_batches(numbers::invalid_unsigned_int) - , n_q_points(numbers::invalid_unsigned_int) - , n_q_points_scalar(numbers::invalid_unsigned_int) - , mapping(&mapping) - , fe(&fe) - , JxW_ptr(nullptr) - , 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) - , is_reinitialized(false) - , is_interior(true) +/** + * This class provides an interface to the evaluation of interpolated solution + * values and gradients on cells on arbitrary reference point positions. These + * points can change from cell to cell, both with respect to their quantity as + * well to the location. The two typical use cases are evaluations on + * non-matching grids and particle simulations. + * + * The use of this class is similar to FEValues or FEEvaluation: The class is + * first initialized to a cell by calling `FEPointEvaluation::reinit(cell, + * unit_points)`, with the main difference to the other concepts that the + * underlying points in reference coordinates need to be passed along. Then, + * upon call to evaluate() or integrate(), the user can compute information at + * the give points. Eventually, the access functions get_value() or + * get_gradient() allow to query this information at a specific point index. + * + * The functionality is similar to creating an FEValues object with a + * Quadrature object on the `unit_points` on every cell separately and then + * calling FEValues::get_function_values or FEValues::get_function_gradients, + * and for some elements and mappings this is what actually happens + * internally. For specific combinations of Mapping and FiniteElement + * realizations, however, there is a much more efficient implementation that + * avoids the memory allocation and other expensive start-up cost of + * FEValues. Currently, the functionality is specialized for mappings derived + * from MappingQ and MappingCartesian and for finite elements with tensor + * product structure that work with the + * @ref matrixfree + * module. In those cases, the cost implied + * by this class is similar (or sometimes even somewhat lower) than using + * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`. + */ +template +class FEPointEvaluation + : public FEPointEvaluationBase { - setup(first_selected_component); -} +public: + static constexpr unsigned int dimension = dim; + static constexpr unsigned int n_components = n_components_; + using number_type = Number; + using ScalarNumber = + typename internal::VectorizedArrayTrait::value_type; + using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; + using ETT = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits; + using value_type = typename ETT::value_type; + using scalar_value_type = typename ETT::scalar_value_type; + using vectorized_value_type = typename ETT::vectorized_value_type; + using gradient_type = typename ETT::gradient_type; + using interface_vectorized_gradient_type = + typename ETT::interface_vectorized_gradient_type; -template -FEPointEvaluationBase:: - FEPointEvaluationBase( + /** + * Constructor. + * + * @param mapping The Mapping class describing the actual geometry of a cell + * passed to the evaluate() function. + * + * @param fe The FiniteElement object that is used for the evaluation, which + * is typically the same on all cells to be evaluated. + * + * @param update_flags Specify the quantities to be computed by the mapping + * during the call of reinit(). During evaluate() or integrate(), this data + * is queried to produce the desired result (e.g., the gradient of a finite + * element solution). + * + * @param first_selected_component For multi-component FiniteElement + * objects, this parameter allows to select a range of `n_components` + * components starting from this parameter. + */ + FEPointEvaluation(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component = 0); + + /** + * Constructor to make the present class able to re-use the geometry + * data also used by other `FEPointEvaluation` objects. + * + * @param mapping_info The MappingInfo class describes the geometry-related + * data for evaluating finite-element solutions. This object enables to + * construct such an object on the outside, possibly re-using it between + * several objects or between several calls to the same cell and unit points. + * + * @param fe The FiniteElement object that is used for the evaluation, which + * is typically the same on all cells to be evaluated. + * + * @param first_selected_component For multi-component FiniteElement + * objects, this parameter allows to select a range of `n_components` + * components starting from this parameter. + */ + FEPointEvaluation( NonMatching::MappingInfo &mapping_info, const FiniteElement &fe, - const unsigned int first_selected_component, - const bool is_interior) - : n_q_batches(numbers::invalid_unsigned_int) - , n_q_points(numbers::invalid_unsigned_int) - , n_q_points_scalar(numbers::invalid_unsigned_int) - , mapping(&mapping_info.get_mapping()) - , fe(&fe) - , JxW_ptr(nullptr) - , 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) - , is_reinitialized(false) - , is_interior(is_interior) -{ - setup(first_selected_component); - connection_is_reinitialized = mapping_info.connect_is_reinitialized( - [this]() { this->is_reinitialized = false; }); -} + const unsigned int first_selected_component = 0); + /** + * Set up the mapping information for the given cell, e.g., by computing the + * Jacobian of the mapping for the given points if gradients of the functions + * are requested. + * + * @param[in] cell An iterator to the current cell + * + * @param[in] unit_points List of points in the reference locations of the + * current cell where the FiniteElement object should be + * evaluated/integrated in the evaluate() and integrate() functions. + */ + void + reinit(const typename Triangulation::cell_iterator &cell, + const ArrayView> &unit_points); + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the single cell in the MappingInfo object. + */ + void + reinit(); -template -FEPointEvaluationBase:: - FEPointEvaluationBase( - FEPointEvaluationBase &other) noexcept - : n_q_batches(other.n_q_batches) - , n_q_points(other.n_q_points) - , n_q_points_scalar(other.n_q_points_scalar) - , mapping(other.mapping) - , fe(other.fe) - , poly(other.poly) - , use_linear_path(other.use_linear_path) - , renumber(other.renumber) - , solution_renumbered(other.solution_renumbered) - , solution_renumbered_vectorized(other.solution_renumbered_vectorized) - , values(other.values) - , gradients(other.gradients) - , dofs_per_component(other.dofs_per_component) - , dofs_per_component_face(other.dofs_per_component_face) - , component_in_base_element(other.component_in_base_element) - , nonzero_shape_function_component(other.nonzero_shape_function_component) - , update_flags(other.update_flags) - , fe_values(other.fe_values) - , mapping_info_on_the_fly( - other.mapping_info_on_the_fly ? - std::make_unique>( - *mapping, - update_flags) : - nullptr) - , mapping_info(other.mapping_info) - , current_cell_index(other.current_cell_index) - , current_face_number(other.current_face_number) - , fast_path(other.fast_path) - , is_reinitialized(false) - , shapes(other.shapes) - , shapes_faces(other.shapes_faces) - , is_interior(other.is_interior) -{ - connection_is_reinitialized = mapping_info->connect_is_reinitialized( - [this]() { this->is_reinitialized = false; }); -} + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the cell in the MappingInfo object. + */ + void + reinit(const unsigned int cell_index); + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element read out by + * `FEEvaluation::read_dof_values(global_vector)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + template + void + evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); -template -FEPointEvaluationBase:: - FEPointEvaluationBase( - FEPointEvaluationBase - &&other) noexcept - : n_q_batches(other.n_q_batches) - , n_q_points(other.n_q_points) - , n_q_points_scalar(other.n_q_points_scalar) - , mapping(other.mapping) - , fe(other.fe) - , poly(other.poly) - , use_linear_path(other.use_linear_path) - , renumber(other.renumber) - , solution_renumbered(other.solution_renumbered) - , solution_renumbered_vectorized(other.solution_renumbered_vectorized) - , values(other.values) - , gradients(other.gradients) - , dofs_per_component(other.dofs_per_component) - , dofs_per_component_face(other.dofs_per_component_face) - , component_in_base_element(other.component_in_base_element) - , nonzero_shape_function_component(other.nonzero_shape_function_component) - , update_flags(other.update_flags) - , fe_values(other.fe_values) - , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly)) - , mapping_info(other.mapping_info) - , current_cell_index(other.current_cell_index) - , current_face_number(other.current_face_number) - , fast_path(other.fast_path) - , is_reinitialized(false) - , shapes(other.shapes) - , shapes_faces(other.shapes_faces) - , is_interior(other.is_interior) -{ - connection_is_reinitialized = mapping_info->connect_is_reinitialized( - [this]() { this->is_reinitialized = false; }); -} - - + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element as returned by `cell->get_dof_values(global_vector, + * solution_values)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + void + evaluate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); -template -FEPointEvaluationBase:: - ~FEPointEvaluationBase() -{ - connection_is_reinitialized.disconnect(); -} + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + template + void + integrate(const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used to during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + void + integrate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points. This is + * similar to the integration of a bilinear form in terms of the test + * function, with the difference that this formula does not include a `JxW` + * factor (in contrast to the integrate function of this class). This allows + * the class to naturally embed point information (e.g. particles) into a + * finite element formulation. + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + template + void + test_and_sum( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); -template -void -FEPointEvaluationBase::setup( - const unsigned int first_selected_component) -{ - AssertIndexRange(first_selected_component + n_components, - fe->n_components() + 1); + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points. This is + * similar to the integration of a bilinear form in terms of the test + * function, with the difference that this formula does not include a `JxW` + * factor (in contrast to the integrate function of this class). This allows + * the class to naturally embed point information (e.g. particles) into a + * finite element formulation. + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + void + test_and_sum(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); - shapes.reserve(100); + /** + * 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, Number> + normal_vector(const unsigned int point_index) const; - bool same_base_element = true; - unsigned int base_element_number = 0; - component_in_base_element = 0; - unsigned int component = 0; - for (; base_element_number < fe->n_base_elements(); ++base_element_number) - if (component + fe->element_multiplicity(base_element_number) > - first_selected_component) - { - if (first_selected_component + n_components > - component + fe->element_multiplicity(base_element_number)) - same_base_element = false; - component_in_base_element = first_selected_component - component; - break; - } - else - component += fe->element_multiplicity(base_element_number); +private: + static constexpr std::size_t n_lanes_user_interface = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t n_lanes_internal = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t stride = + internal::VectorizedArrayTrait::stride(); - if (internal::FEPointEvaluation::is_fast_path_supported(*mapping) && - internal::FEPointEvaluation::is_fast_path_supported( - *fe, base_element_number) && - same_base_element) - { - shape_info.reinit(QMidpoint<1>(), *fe, base_element_number); - renumber = shape_info.lexicographic_numbering; - dofs_per_component = shape_info.dofs_per_component_on_cell; - dofs_per_component_face = shape_info.dofs_per_component_on_face; - poly = internal::FEPointEvaluation::get_polynomial_space( - fe->base_element(base_element_number)); + /** + * Resizes necessary data fields, reads in and renumbers solution values. + * Interpolates onto face if face path is selected. + */ + template + void + prepare_evaluate_fast( + const StridedArrayView &solution_values); - bool is_lexicographic = true; - for (unsigned int i = 0; i < renumber.size(); ++i) - if (i != renumber[i]) - is_lexicographic = false; + /** + * Evaluates the actual interpolation on the cell or face for a quadrature + * batch. + */ + template + void + compute_evaluate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags, + const unsigned int n_shapes, + const unsigned int qb, + vectorized_value_type &value, + interface_vectorized_gradient_type &gradient); - if (is_lexicographic) - renumber.clear(); + /** + * Fast path of the evaluate function. + */ + template + void + evaluate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); - use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. && - poly[0].value(1.) == 0. && poly[1].value(0.) == 0. && - poly[1].value(1.) == 1.) && - (fe->n_components() == n_components); + /** + * Slow path of the evaluate function using FEValues. + */ + template + void + evaluate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); - const unsigned int size_face = 3 * dofs_per_component_face * n_components; - const unsigned int size_cell = dofs_per_component * n_components; - scratch_data_scalar.resize(size_face + size_cell); + /** + * Integrates the product of the data passed in by submit_value() and + * submit_gradient() with the values or gradients of test functions on the + * cell or face for a given quadrature batch. + */ + template + void + compute_integrate_fast( + const EvaluationFlags::EvaluationFlags &integration_flags, + const unsigned int n_shapes, + const unsigned int qb, + const vectorized_value_type value, + const interface_vectorized_gradient_type gradient, + vectorized_value_type *solution_values_vectorized_linear); - solution_renumbered.resize(dofs_per_component); - solution_renumbered_vectorized.resize(dofs_per_component); + /** + * Addition across the lanes of VectorizedArray as accumulated by the + * compute_integrate_fast_function(), writing the sum into the result vector. + * Applies face contributions to cell contributions for face path. + */ + template + void + finish_integrate_fast( + const StridedArrayView &solution_values, + vectorized_value_type *solution_values_vectorized_linear, + const bool sum_into_values); - fast_path = true; - } - else - { - nonzero_shape_function_component.resize(fe->n_dofs_per_cell()); - for (unsigned int d = 0; d < n_components; ++d) - { - const unsigned int component = first_selected_component + d; - for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) - { - const bool is_primitive = - fe->is_primitive() || fe->is_primitive(i); - if (is_primitive) - nonzero_shape_function_component[i][d] = - (component == fe->system_to_component_index(i).first); - else - nonzero_shape_function_component[i][d] = - (fe->get_nonzero_components(i)[component] == true); - } - } + /** + * Fast path of the integrate function. + */ + template + void + integrate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); - fast_path = false; - } -} + /** + * Slow path of the integrate function using FEValues. + */ + template + void + integrate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); + /** + * Implementation of the integrate/test_and_sum function. + */ + template + void + do_integrate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); +}; -template -template -inline void -FEPointEvaluationBase::do_reinit() -{ - const unsigned int geometry_index = - mapping_info->template compute_geometry_index_offset( - current_cell_index, current_face_number); - cell_type = mapping_info->get_cell_type(geometry_index); +/** + * This class provides an interface to the evaluation of interpolated solution + * values and gradients on faces on arbitrary reference point positions. These + * points can change from face to face, both with respect to their quantity as + * well to the location. A typical use case is evaluations on non-matching + * grids. + * + * The use of this class is similar to FEEvaluation: In the constructor, a + * reference to a NonMatching::MappingInfo object is passed, where the + * quadrature points in reference position is stored together with the mapping + * information. The class is then reinitialized to a cell by calling + * `FEFacePointEvaluation::reinit(face_index)` or + * `FEFacePointEvaluation::reinit(cell_index, face_number)`. Then, upon call to + * evaluate() or integrate(), the user can compute information at the given + * points. Eventually, the access functions get_value() or get_gradient() allow + * to query this information at a specific point index. + */ +template +class FEFacePointEvaluation + : public FEPointEvaluationBase +{ +public: + static constexpr unsigned int dimension = dim; + static constexpr unsigned int n_components = n_components_; - const_cast(n_q_points_scalar) = - mapping_info->get_n_q_points_unvectorized(geometry_index); + using number_type = Number; - // round up n_q_points_scalar / n_lanes_internal - const_cast(n_q_batches) = - (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal; + using ScalarNumber = + typename internal::VectorizedArrayTrait::value_type; + using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; + using ETT = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits; + using value_type = typename ETT::value_type; + using scalar_value_type = typename ETT::scalar_value_type; + using vectorized_value_type = typename ETT::vectorized_value_type; + using gradient_type = typename ETT::gradient_type; + using interface_vectorized_gradient_type = + typename ETT::interface_vectorized_gradient_type; - const unsigned int n_q_points_before = n_q_points; + /** + * Constructor. Allows to select if interior or exterior face is selected. + */ + FEFacePointEvaluation( + NonMatching::MappingInfo &mapping_info, + const FiniteElement &fe, + const bool is_interior = true, + const unsigned int first_selected_component = 0); - const_cast(n_q_points) = - (stride == 1) ? n_q_batches : n_q_points_scalar; + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the face in the MappingInfo object. Used in element-centric loops (ECL). + */ + void + reinit(const unsigned int cell_index, const unsigned int face_number); - if (n_q_points != n_q_points_before) - { - if (update_flags & update_values) - values.resize(n_q_points); - if (update_flags & update_gradients) - gradients.resize(n_q_points); - } + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the face in the MappingInfo object. Used in face-centric loops (FCL). + */ + void + reinit(const unsigned int face_index); - if (n_q_points == 0) - { - is_reinitialized = true; - return; - } + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element read out by + * `FEEvaluation::read_dof_values(global_vector)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + template + void + evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); - // set unit point pointer - const unsigned int unit_point_offset = - mapping_info->compute_unit_point_index_offset(geometry_index); + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element as returned by `cell->get_dof_values(global_vector, + * solution_values)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + void + evaluate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); - if (is_face) - unit_point_faces_ptr = - mapping_info->get_unit_point_faces(unit_point_offset); - else - unit_point_ptr = mapping_info->get_unit_point(unit_point_offset); + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + template + void + integrate(const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); - // set data pointers - const unsigned int data_offset = - mapping_info->compute_data_index_offset(geometry_index); - const unsigned int compressed_data_offset = - mapping_info->compute_compressed_data_index_offset(geometry_index); -#ifdef DEBUG - const UpdateFlags update_flags_mapping = - mapping_info->get_update_flags_mapping(); - if (update_flags_mapping & UpdateFlags::update_quadrature_points) - real_point_ptr = mapping_info->get_real_point(data_offset); - if (update_flags_mapping & UpdateFlags::update_jacobians) - jacobian_ptr = - mapping_info->get_jacobian(compressed_data_offset, is_interior); - if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) - inverse_jacobian_ptr = - mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior); - if (update_flags_mapping & UpdateFlags::update_normal_vectors) - normal_ptr = mapping_info->get_normal_vector(data_offset); - if (update_flags_mapping & UpdateFlags::update_JxW_values) - JxW_ptr = mapping_info->get_JxW(data_offset); -#else - real_point_ptr = mapping_info->get_real_point(data_offset); - jacobian_ptr = - mapping_info->get_jacobian(compressed_data_offset, is_interior); - inverse_jacobian_ptr = - mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior); - normal_ptr = mapping_info->get_normal_vector(data_offset); - JxW_ptr = mapping_info->get_JxW(data_offset); -#endif + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used to during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + void + integrate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); - if (!is_linear && fast_path) - { - const std::size_t n_shapes = poly.size(); - - for (unsigned int qb = 0; qb < n_q_batches; ++qb) - if (is_face) - { - if (dim > 1) - { - shapes_faces.resize_fast(n_q_batches * n_shapes); - internal::compute_values_of_array( - shapes_faces.data() + qb * n_shapes, - poly, - unit_point_faces_ptr[qb], - update_flags & UpdateFlags::update_gradients ? 1 : 0); - } - } - else - { - shapes.resize_fast(n_q_batches * n_shapes); - internal::compute_values_of_array( - shapes.data() + qb * n_shapes, - poly, - unit_point_ptr[qb], - update_flags & UpdateFlags::update_gradients ? 1 : 0); - } - } - - is_reinitialized = true; -} - - - -template -inline const typename FEPointEvaluationBase::value_type & -FEPointEvaluationBase::get_value( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, values.size()); - return values[point_index]; -} + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + template + void + test_and_sum( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used to during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. For the default value + * `sum_into_values=false` every value of @p solution_values is zeroed out. + * + */ + void + test_and_sum(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + /** + * Evaluate values and gradients in face for the selected face (lane) of the + * batch. Default stride into the face dofs is width of + * VectorizedArray which is the default + * vectorization over faces for FEFaceEvaluation. + */ + template + void + evaluate_in_face(const ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); -template -inline const typename FEPointEvaluationBase::gradient_type & -FEPointEvaluationBase::get_gradient( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, gradients.size()); - return gradients[point_index]; -} + /** + * Integrate values and gradients in face for the selected face (lane) of the + * batch. Default stride into the face dofs is width of + * VectorizedArray which is the default + * vectorization over faces for FEFaceEvaluation. + */ + template + void + integrate_in_face(ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + /** + * 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, Number> + normal_vector(const unsigned int point_index) const; +private: + static constexpr std::size_t n_lanes_user_interface = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t n_lanes_internal = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t stride = + internal::VectorizedArrayTrait::stride(); -template -inline void -FEPointEvaluationBase::submit_value( - const value_type &value, - const unsigned int point_index) -{ - AssertIndexRange(point_index, n_q_points); - values[point_index] = value; -} + template + void + do_evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + template + void + do_integrate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); + /** + * Actually does the evaluation templated on the chosen code path (linear or + * higher order). + */ + template + void + do_evaluate_in_face(const ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); -template -inline void -FEPointEvaluationBase::submit_gradient( - const gradient_type &gradient, - const unsigned int point_index) -{ - AssertIndexRange(point_index, n_q_points); - gradients[point_index] = gradient; -} + /** + * Actually does the integration templated on the chosen code path (linear or + * higher order). + */ + template + void + do_integrate_in_face( + ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); +}; +// ----------------------- template and inline function ---------------------- template -inline DerivativeForm<1, dim, spacedim, Number> -FEPointEvaluationBase::jacobian( - const unsigned int point_index) const +FEPointEvaluationBase:: + FEPointEvaluationBase(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component) + : n_q_batches(numbers::invalid_unsigned_int) + , n_q_points(numbers::invalid_unsigned_int) + , n_q_points_scalar(numbers::invalid_unsigned_int) + , mapping(&mapping) + , fe(&fe) + , JxW_ptr(nullptr) + , 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) + , is_reinitialized(false) + , is_interior(true) { - AssertIndexRange(point_index, n_q_points); - Assert(jacobian_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_jacobians")); - return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions:: - GeometryType::affine ? - 0 : - point_index]; + setup(first_selected_component); } template -inline DerivativeForm<1, spacedim, dim, Number> -FEPointEvaluationBase::inverse_jacobian( - const unsigned int point_index) const +FEPointEvaluationBase:: + FEPointEvaluationBase( + NonMatching::MappingInfo &mapping_info, + const FiniteElement &fe, + const unsigned int first_selected_component, + const bool is_interior) + : n_q_batches(numbers::invalid_unsigned_int) + , n_q_points(numbers::invalid_unsigned_int) + , n_q_points_scalar(numbers::invalid_unsigned_int) + , mapping(&mapping_info.get_mapping()) + , fe(&fe) + , JxW_ptr(nullptr) + , 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) + , is_reinitialized(false) + , is_interior(is_interior) { - AssertIndexRange(point_index, n_q_points); - Assert(inverse_jacobian_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_inverse_jacobians")); - return inverse_jacobian_ptr - [cell_type <= - ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ? - 0 : - point_index]; + setup(first_selected_component); + connection_is_reinitialized = mapping_info.connect_is_reinitialized( + [this]() { this->is_reinitialized = false; }); } template -inline Number -FEPointEvaluationBase::JxW( - const unsigned int point_index) const +FEPointEvaluationBase:: + FEPointEvaluationBase( + FEPointEvaluationBase &other) noexcept + : n_q_batches(other.n_q_batches) + , n_q_points(other.n_q_points) + , n_q_points_scalar(other.n_q_points_scalar) + , mapping(other.mapping) + , fe(other.fe) + , poly(other.poly) + , use_linear_path(other.use_linear_path) + , renumber(other.renumber) + , solution_renumbered(other.solution_renumbered) + , solution_renumbered_vectorized(other.solution_renumbered_vectorized) + , values(other.values) + , gradients(other.gradients) + , dofs_per_component(other.dofs_per_component) + , dofs_per_component_face(other.dofs_per_component_face) + , component_in_base_element(other.component_in_base_element) + , nonzero_shape_function_component(other.nonzero_shape_function_component) + , update_flags(other.update_flags) + , fe_values(other.fe_values) + , mapping_info_on_the_fly( + other.mapping_info_on_the_fly ? + std::make_unique>( + *mapping, + update_flags) : + nullptr) + , mapping_info(other.mapping_info) + , current_cell_index(other.current_cell_index) + , current_face_number(other.current_face_number) + , fast_path(other.fast_path) + , is_reinitialized(false) + , shapes(other.shapes) + , shapes_faces(other.shapes_faces) + , is_interior(other.is_interior) { - AssertIndexRange(point_index, n_q_points); - Assert(JxW_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_JxW_values")); - return JxW_ptr[point_index]; + connection_is_reinitialized = mapping_info->connect_is_reinitialized( + [this]() { this->is_reinitialized = false; }); } template -inline Point -FEPointEvaluationBase::real_point( - const unsigned int point_index) const +FEPointEvaluationBase:: + FEPointEvaluationBase( + FEPointEvaluationBase + &&other) noexcept + : n_q_batches(other.n_q_batches) + , n_q_points(other.n_q_points) + , n_q_points_scalar(other.n_q_points_scalar) + , mapping(other.mapping) + , fe(other.fe) + , poly(other.poly) + , use_linear_path(other.use_linear_path) + , renumber(other.renumber) + , solution_renumbered(other.solution_renumbered) + , solution_renumbered_vectorized(other.solution_renumbered_vectorized) + , values(other.values) + , gradients(other.gradients) + , dofs_per_component(other.dofs_per_component) + , dofs_per_component_face(other.dofs_per_component_face) + , component_in_base_element(other.component_in_base_element) + , nonzero_shape_function_component(other.nonzero_shape_function_component) + , update_flags(other.update_flags) + , fe_values(other.fe_values) + , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly)) + , mapping_info(other.mapping_info) + , current_cell_index(other.current_cell_index) + , current_face_number(other.current_face_number) + , fast_path(other.fast_path) + , is_reinitialized(false) + , shapes(other.shapes) + , shapes_faces(other.shapes_faces) + , is_interior(other.is_interior) { - AssertIndexRange(point_index, n_q_points); - Assert(real_point_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_quadrature_points")); - return real_point_ptr[point_index]; + connection_is_reinitialized = mapping_info->connect_is_reinitialized( + [this]() { this->is_reinitialized = false; }); } template -inline Point -FEPointEvaluationBase::unit_point( - const unsigned int point_index) const +FEPointEvaluationBase:: + ~FEPointEvaluationBase() { - AssertIndexRange(point_index, n_q_points); - Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!")); - Point unit_point; - for (unsigned int d = 0; d < dim; ++d) - unit_point[d] = internal::VectorizedArrayTrait::get_from_vectorized( - unit_point_ptr[point_index / stride][d], point_index % stride); - return unit_point; + connection_is_reinitialized.disconnect(); } template -inline std_cxx20::ranges::iota_view -FEPointEvaluationBase:: - quadrature_point_indices() const +void +FEPointEvaluationBase::setup( + const unsigned int first_selected_component) { - return {0U, n_q_points}; -} + AssertIndexRange(first_selected_component + n_components, + fe->n_components() + 1); + shapes.reserve(100); + bool same_base_element = true; + unsigned int base_element_number = 0; + component_in_base_element = 0; + unsigned int component = 0; + for (; base_element_number < fe->n_base_elements(); ++base_element_number) + if (component + fe->element_multiplicity(base_element_number) > + first_selected_component) + { + if (first_selected_component + n_components > + component + fe->element_multiplicity(base_element_number)) + same_base_element = false; + component_in_base_element = first_selected_component - component; + break; + } + else + component += fe->element_multiplicity(base_element_number); -/** - * This class provides an interface to the evaluation of interpolated solution - * values and gradients on cells on arbitrary reference point positions. These - * points can change from cell to cell, both with respect to their quantity as - * well to the location. The two typical use cases are evaluations on - * non-matching grids and particle simulations. - * - * The use of this class is similar to FEValues or FEEvaluation: The class is - * first initialized to a cell by calling `FEPointEvaluation::reinit(cell, - * unit_points)`, with the main difference to the other concepts that the - * underlying points in reference coordinates need to be passed along. Then, - * upon call to evaluate() or integrate(), the user can compute information at - * the give points. Eventually, the access functions get_value() or - * get_gradient() allow to query this information at a specific point index. - * - * The functionality is similar to creating an FEValues object with a - * Quadrature object on the `unit_points` on every cell separately and then - * calling FEValues::get_function_values or FEValues::get_function_gradients, - * and for some elements and mappings this is what actually happens - * internally. For specific combinations of Mapping and FiniteElement - * realizations, however, there is a much more efficient implementation that - * avoids the memory allocation and other expensive start-up cost of - * FEValues. Currently, the functionality is specialized for mappings derived - * from MappingQ and MappingCartesian and for finite elements with tensor - * product structure that work with the - * @ref matrixfree - * module. In those cases, the cost implied - * by this class is similar (or sometimes even somewhat lower) than using - * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`. - */ -template -class FEPointEvaluation - : public FEPointEvaluationBase -{ -public: - static constexpr unsigned int dimension = dim; - static constexpr unsigned int n_components = n_components_; + if (internal::FEPointEvaluation::is_fast_path_supported(*mapping) && + internal::FEPointEvaluation::is_fast_path_supported( + *fe, base_element_number) && + same_base_element) + { + shape_info.reinit(QMidpoint<1>(), *fe, base_element_number); + renumber = shape_info.lexicographic_numbering; + dofs_per_component = shape_info.dofs_per_component_on_cell; + dofs_per_component_face = shape_info.dofs_per_component_on_face; + poly = internal::FEPointEvaluation::get_polynomial_space( + fe->base_element(base_element_number)); - using number_type = Number; + bool is_lexicographic = true; + for (unsigned int i = 0; i < renumber.size(); ++i) + if (i != renumber[i]) + is_lexicographic = false; - using ScalarNumber = - typename internal::VectorizedArrayTrait::value_type; - using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< - Number>::vectorized_value_type; - using ETT = typename internal::FEPointEvaluation:: - EvaluatorTypeTraits; - using value_type = typename ETT::value_type; - using scalar_value_type = typename ETT::scalar_value_type; - using vectorized_value_type = typename ETT::vectorized_value_type; - using gradient_type = typename ETT::gradient_type; - using interface_vectorized_gradient_type = - typename ETT::interface_vectorized_gradient_type; + if (is_lexicographic) + renumber.clear(); - /** - * Constructor. - * - * @param mapping The Mapping class describing the actual geometry of a cell - * passed to the evaluate() function. - * - * @param fe The FiniteElement object that is used for the evaluation, which - * is typically the same on all cells to be evaluated. - * - * @param update_flags Specify the quantities to be computed by the mapping - * during the call of reinit(). During evaluate() or integrate(), this data - * is queried to produce the desired result (e.g., the gradient of a finite - * element solution). - * - * @param first_selected_component For multi-component FiniteElement - * objects, this parameter allows to select a range of `n_components` - * components starting from this parameter. - */ - FEPointEvaluation(const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component = 0) - : FEPointEvaluationBase( - mapping, - fe, - update_flags, - first_selected_component) - {} + use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. && + poly[0].value(1.) == 0. && poly[1].value(0.) == 0. && + poly[1].value(1.) == 1.) && + (fe->n_components() == n_components); - /** - * Constructor to make the present class able to re-use the geometry - * data also used by other `FEPointEvaluation` objects. - * - * @param mapping_info The MappingInfo class describes the geometry-related - * data for evaluating finite-element solutions. This object enables to - * construct such an object on the outside, possibly re-using it between - * several objects or between several calls to the same cell and unit points. - * - * @param fe The FiniteElement object that is used for the evaluation, which - * is typically the same on all cells to be evaluated. - * - * @param first_selected_component For multi-component FiniteElement - * objects, this parameter allows to select a range of `n_components` - * components starting from this parameter. - */ - FEPointEvaluation( - NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, - const unsigned int first_selected_component = 0) - : FEPointEvaluationBase( - mapping_info, - fe, - first_selected_component) - {} + const unsigned int size_face = 3 * dofs_per_component_face * n_components; + const unsigned int size_cell = dofs_per_component * n_components; + scratch_data_scalar.resize(size_face + size_cell); - /** - * Set up the mapping information for the given cell, e.g., by computing the - * Jacobian of the mapping for the given points if gradients of the functions - * are requested. - * - * @param[in] cell An iterator to the current cell - * - * @param[in] unit_points List of points in the reference locations of the - * current cell where the FiniteElement object should be - * evaluated/integrated in the evaluate() and integrate() functions. - */ - void - reinit(const typename Triangulation::cell_iterator &cell, - const ArrayView> &unit_points); - - /** - * Reinitialize the evaluator to point to the correct precomputed mapping of - * the single cell in the MappingInfo object. - */ - void - reinit(); + solution_renumbered.resize(dofs_per_component); + solution_renumbered_vectorized.resize(dofs_per_component); - /** - * Reinitialize the evaluator to point to the correct precomputed mapping of - * the cell in the MappingInfo object. - */ - void - reinit(const unsigned int cell_index); + fast_path = true; + } + else + { + nonzero_shape_function_component.resize(fe->n_dofs_per_cell()); + for (unsigned int d = 0; d < n_components; ++d) + { + const unsigned int component = first_selected_component + d; + for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) + { + const bool is_primitive = + fe->is_primitive() || fe->is_primitive(i); + if (is_primitive) + nonzero_shape_function_component[i][d] = + (component == fe->system_to_component_index(i).first); + else + nonzero_shape_function_component[i][d] = + (fe->get_nonzero_components(i)[component] == true); + } + } + fast_path = false; + } +} - /** - * This function interpolates the finite element solution, represented by - * `solution_values`, on the cell and `unit_points` passed to reinit(). - * - * @param[in] solution_values This array is supposed to contain the unknown - * values on the element read out by - * `FEEvaluation::read_dof_values(global_vector)`. - * - * @param[in] evaluation_flags Flags specifying which quantities should be - * evaluated at the points. - */ - template - void - evaluate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - /** - * This function interpolates the finite element solution, represented by - * `solution_values`, on the cell and `unit_points` passed to reinit(). - * - * @param[in] solution_values This array is supposed to contain the unknown - * values on the element as returned by `cell->get_dof_values(global_vector, - * solution_values)`. - * - * @param[in] evaluation_flags Flags specifying which quantities should be - * evaluated at the points. - */ - void - evaluate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `FEEvaluation::set_dof_values(global_vector)` or - * `FEEvaluation::distribute_local_to_global(global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - template - void - integrate(const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); +template +template +inline void +FEPointEvaluationBase::do_reinit() +{ + const unsigned int geometry_index = + mapping_info->template compute_geometry_index_offset( + current_cell_index, current_face_number); - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used to during - * `cell->set_dof_values(solution_values, global_vector)` or - * `cell->distribute_local_to_global(solution_values, global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - void - integrate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); + cell_type = mapping_info->get_cell_type(geometry_index); - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points. This is - * similar to the integration of a bilinear form in terms of the test - * function, with the difference that this formula does not include a `JxW` - * factor (in contrast to the integrate function of this class). This allows - * the class to naturally embed point information (e.g. particles) into a - * finite element formulation. - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `FEEvaluation::set_dof_values(global_vector)` or - * `FEEvaluation::distribute_local_to_global(global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - template - void - test_and_sum( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); + const_cast(n_q_points_scalar) = + mapping_info->get_n_q_points_unvectorized(geometry_index); - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points. This is - * similar to the integration of a bilinear form in terms of the test - * function, with the difference that this formula does not include a `JxW` - * factor (in contrast to the integrate function of this class). This allows - * the class to naturally embed point information (e.g. particles) into a - * finite element formulation. - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `cell->set_dof_values(solution_values, global_vector)` or - * `cell->distribute_local_to_global(solution_values, global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - void - test_and_sum(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); + // round up n_q_points_scalar / n_lanes_internal + const_cast(n_q_batches) = + (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal; - /** - * 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, Number> - normal_vector(const unsigned int point_index) const; + const unsigned int n_q_points_before = n_q_points; -private: - static constexpr std::size_t n_lanes_user_interface = - internal::VectorizedArrayTrait::width(); - static constexpr std::size_t n_lanes_internal = - internal::VectorizedArrayTrait::width(); - static constexpr std::size_t stride = - internal::VectorizedArrayTrait::stride(); + const_cast(n_q_points) = + (stride == 1) ? n_q_batches : n_q_points_scalar; - /** - * Resizes necessary data fields, reads in and renumbers solution values. - * Interpolates onto face if face path is selected. - */ - template - void - prepare_evaluate_fast( - const StridedArrayView &solution_values); + if (n_q_points != n_q_points_before) + { + if (update_flags & update_values) + values.resize(n_q_points); + if (update_flags & update_gradients) + gradients.resize(n_q_points); + } - /** - * Evaluates the actual interpolation on the cell or face for a quadrature - * batch. - */ - template - void - compute_evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags, - const unsigned int n_shapes, - const unsigned int qb, - vectorized_value_type &value, - interface_vectorized_gradient_type &gradient); + if (n_q_points == 0) + { + is_reinitialized = true; + return; + } - /** - * Fast path of the evaluate function. - */ - template - void - evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); + // set unit point pointer + const unsigned int unit_point_offset = + mapping_info->compute_unit_point_index_offset(geometry_index); - /** - * Slow path of the evaluate function using FEValues. - */ - template - void - evaluate_slow( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * Integrates the product of the data passed in by submit_value() and - * submit_gradient() with the values or gradients of test functions on the - * cell or face for a given quadrature batch. - */ - template - void - compute_integrate_fast( - const EvaluationFlags::EvaluationFlags &integration_flags, - const unsigned int n_shapes, - const unsigned int qb, - const vectorized_value_type value, - const interface_vectorized_gradient_type gradient, - vectorized_value_type *solution_values_vectorized_linear); + if (is_face) + unit_point_faces_ptr = + mapping_info->get_unit_point_faces(unit_point_offset); + else + unit_point_ptr = mapping_info->get_unit_point(unit_point_offset); - /** - * Addition across the lanes of VectorizedArray as accumulated by the - * compute_integrate_fast_function(), writing the sum into the result vector. - * Applies face contributions to cell contributions for face path. - */ - template - void - finish_integrate_fast( - const StridedArrayView &solution_values, - vectorized_value_type *solution_values_vectorized_linear, - const bool sum_into_values); + // set data pointers + const unsigned int data_offset = + mapping_info->compute_data_index_offset(geometry_index); + const unsigned int compressed_data_offset = + mapping_info->compute_compressed_data_index_offset(geometry_index); +#ifdef DEBUG + const UpdateFlags update_flags_mapping = + mapping_info->get_update_flags_mapping(); + if (update_flags_mapping & UpdateFlags::update_quadrature_points) + real_point_ptr = mapping_info->get_real_point(data_offset); + if (update_flags_mapping & UpdateFlags::update_jacobians) + jacobian_ptr = + mapping_info->get_jacobian(compressed_data_offset, is_interior); + if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) + inverse_jacobian_ptr = + mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior); + if (update_flags_mapping & UpdateFlags::update_normal_vectors) + normal_ptr = mapping_info->get_normal_vector(data_offset); + if (update_flags_mapping & UpdateFlags::update_JxW_values) + JxW_ptr = mapping_info->get_JxW(data_offset); +#else + real_point_ptr = mapping_info->get_real_point(data_offset); + jacobian_ptr = + mapping_info->get_jacobian(compressed_data_offset, is_interior); + inverse_jacobian_ptr = + mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior); + normal_ptr = mapping_info->get_normal_vector(data_offset); + JxW_ptr = mapping_info->get_JxW(data_offset); +#endif - /** - * Fast path of the integrate function. - */ - template - void - integrate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); + if (!is_linear && fast_path) + { + const std::size_t n_shapes = poly.size(); - /** - * Slow path of the integrate function using FEValues. - */ - template - void - integrate_slow( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); + for (unsigned int qb = 0; qb < n_q_batches; ++qb) + if (is_face) + { + if (dim > 1) + { + shapes_faces.resize_fast(n_q_batches * n_shapes); + internal::compute_values_of_array( + shapes_faces.data() + qb * n_shapes, + poly, + unit_point_faces_ptr[qb], + update_flags & UpdateFlags::update_gradients ? 1 : 0); + } + } + else + { + shapes.resize_fast(n_q_batches * n_shapes); + internal::compute_values_of_array( + shapes.data() + qb * n_shapes, + poly, + unit_point_ptr[qb], + update_flags & UpdateFlags::update_gradients ? 1 : 0); + } + } - /** - * Implementation of the integrate/test_and_sum function. - */ - template - void - do_integrate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); -}; + is_reinitialized = true; +} template -inline void -FEPointEvaluation::reinit() +inline const typename FEPointEvaluationBase::value_type & +FEPointEvaluationBase::get_value( + const unsigned int point_index) const { - this->current_cell_index = numbers::invalid_unsigned_int; - this->current_face_number = numbers::invalid_unsigned_int; - - if (this->use_linear_path) - this->template do_reinit(); - else - this->template do_reinit(); + AssertIndexRange(point_index, values.size()); + return values[point_index]; } template -inline void -FEPointEvaluation::reinit( - const typename Triangulation::cell_iterator &cell, - const ArrayView> &unit_points) +inline const typename FEPointEvaluationBase::gradient_type & +FEPointEvaluationBase::get_gradient( + const unsigned int point_index) const { - // reinit is only allowed for mapping computation on the fly - AssertThrow(this->mapping_info_on_the_fly.get() != nullptr, - ExcNotImplemented()); + AssertIndexRange(point_index, gradients.size()); + return gradients[point_index]; +} - this->mapping_info->reinit(cell, unit_points); - if (!this->fast_path) - { - this->fe_values = std::make_shared>( - *this->mapping, - *this->fe, - Quadrature( - std::vector>(unit_points.begin(), unit_points.end())), - this->update_flags); - this->fe_values->reinit(cell); - } - if (this->use_linear_path) - this->template do_reinit(); - else - this->template do_reinit(); +template +inline void +FEPointEvaluationBase::submit_value( + const value_type &value, + const unsigned int point_index) +{ + AssertIndexRange(point_index, n_q_points); + values[point_index] = value; } template inline void -FEPointEvaluation::reinit( - const unsigned int cell_index) +FEPointEvaluationBase::submit_gradient( + const gradient_type &gradient, + const unsigned int point_index) { - this->current_cell_index = cell_index; - this->current_face_number = numbers::invalid_unsigned_int; - - if (this->use_linear_path) - this->template do_reinit(); - else - this->template do_reinit(); - - if (!this->fast_path) - { - std::vector> unit_points(this->n_q_points_scalar); + AssertIndexRange(point_index, n_q_points); + gradients[point_index] = gradient; +} - for (unsigned int v = 0; v < this->n_q_points_scalar; ++v) - for (unsigned int d = 0; d < dim; ++d) - unit_points[v][d] = - this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal]; - this->fe_values = std::make_shared>( - *this->mapping, - *this->fe, - Quadrature( - std::vector>(unit_points.begin(), unit_points.end())), - this->update_flags); - this->fe_values->reinit( - this->mapping_info->get_cell_iterator(this->current_cell_index)); - } +template +inline DerivativeForm<1, dim, spacedim, Number> +FEPointEvaluationBase::jacobian( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(jacobian_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_jacobians")); + return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions:: + GeometryType::affine ? + 0 : + point_index]; } template -template -void -FEPointEvaluation::evaluate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) +inline DerivativeForm<1, spacedim, dim, Number> +FEPointEvaluationBase::inverse_jacobian( + const unsigned int point_index) const { - if (!this->is_reinitialized) - reinit(); - - if (this->n_q_points == 0) - return; + AssertIndexRange(point_index, n_q_points); + Assert(inverse_jacobian_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_inverse_jacobians")); + return inverse_jacobian_ptr + [cell_type <= + ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ? + 0 : + point_index]; +} - Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented()); - if (!((evaluation_flags & EvaluationFlags::values) || - (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags - return; - AssertDimension(solution_values.size(), this->fe->dofs_per_cell); - if (this->fast_path) +template +inline Number +FEPointEvaluationBase::JxW( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(JxW_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_JxW_values")); + return JxW_ptr[point_index]; +} + + + +template +inline Point +FEPointEvaluationBase::real_point( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(real_point_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_quadrature_points")); + return real_point_ptr[point_index]; +} + + + +template +inline Point +FEPointEvaluationBase::unit_point( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!")); + Point unit_point; + for (unsigned int d = 0; d < dim; ++d) + unit_point[d] = internal::VectorizedArrayTrait::get_from_vectorized( + unit_point_ptr[point_index / stride][d], point_index % stride); + return unit_point; +} + + + +template +inline std_cxx20::ranges::iota_view +FEPointEvaluationBase:: + quadrature_point_indices() const +{ + return {0U, n_q_points}; +} + + + +template +FEPointEvaluation::FEPointEvaluation( + NonMatching::MappingInfo &mapping_info, + const FiniteElement &fe, + const unsigned int first_selected_component) + : FEPointEvaluationBase( + mapping_info, + fe, + first_selected_component) +{} + + + +template +FEPointEvaluation::FEPointEvaluation( + const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component) + : FEPointEvaluationBase( + mapping, + fe, + update_flags, + first_selected_component) +{} + + + +template +inline void +FEPointEvaluation::reinit() +{ + this->current_cell_index = numbers::invalid_unsigned_int; + this->current_face_number = numbers::invalid_unsigned_int; + + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); +} + + + +template +inline void +FEPointEvaluation::reinit( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &unit_points) +{ + // reinit is only allowed for mapping computation on the fly + AssertThrow(this->mapping_info_on_the_fly.get() != nullptr, + ExcNotImplemented()); + + this->mapping_info->reinit(cell, unit_points); + + if (!this->fast_path) + { + this->fe_values = std::make_shared>( + *this->mapping, + *this->fe, + Quadrature( + std::vector>(unit_points.begin(), unit_points.end())), + this->update_flags); + this->fe_values->reinit(cell); + } + + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); +} + + + +template +inline void +FEPointEvaluation::reinit( + const unsigned int cell_index) +{ + this->current_cell_index = cell_index; + this->current_face_number = numbers::invalid_unsigned_int; + + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); + + if (!this->fast_path) + { + std::vector> unit_points(this->n_q_points_scalar); + + for (unsigned int v = 0; v < this->n_q_points_scalar; ++v) + for (unsigned int d = 0; d < dim; ++d) + unit_points[v][d] = + this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal]; + + this->fe_values = std::make_shared>( + *this->mapping, + *this->fe, + Quadrature( + std::vector>(unit_points.begin(), unit_points.end())), + this->update_flags); + + this->fe_values->reinit( + this->mapping_info->get_cell_iterator(this->current_cell_index)); + } +} + + + +template +template +void +FEPointEvaluation::evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) +{ + if (!this->is_reinitialized) + reinit(); + + if (this->n_q_points == 0) + return; + + Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented()); + + if (!((evaluation_flags & EvaluationFlags::values) || + (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags + return; + + AssertDimension(solution_values.size(), this->fe->dofs_per_cell); + if (this->fast_path) { if (this->use_linear_path) evaluate_fast(solution_values, evaluation_flags); @@ -2766,286 +3071,6 @@ FEPointEvaluation::normal_vector( -/** - * This class provides an interface to the evaluation of interpolated solution - * values and gradients on faces on arbitrary reference point positions. These - * points can change from face to face, both with respect to their quantity as - * well to the location. A typical use case is evaluations on non-matching - * grids. - * - * The use of this class is similar to FEEvaluation: In the constructor, a - * reference to a NonMatching::MappingInfo object is passed, where the - * quadrature points in reference position is stored together with the mapping - * information. The class is then reinitialized to a cell by calling - * `FEFacePointEvaluation::reinit(face_index)` or - * `FEFacePointEvaluation::reinit(cell_index, face_number)`. Then, upon call to - * evaluate() or integrate(), the user can compute information at the given - * points. Eventually, the access functions get_value() or get_gradient() allow - * to query this information at a specific point index. - */ -template -class FEFacePointEvaluation - : public FEPointEvaluationBase -{ -public: - static constexpr unsigned int dimension = dim; - static constexpr unsigned int n_components = n_components_; - - using number_type = Number; - - using ScalarNumber = - typename internal::VectorizedArrayTrait::value_type; - using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< - Number>::vectorized_value_type; - using ETT = typename internal::FEPointEvaluation:: - EvaluatorTypeTraits; - using value_type = typename ETT::value_type; - using scalar_value_type = typename ETT::scalar_value_type; - using vectorized_value_type = typename ETT::vectorized_value_type; - using gradient_type = typename ETT::gradient_type; - using interface_vectorized_gradient_type = - typename ETT::interface_vectorized_gradient_type; - - /** - * Constructor. Allows to select if interior or exterior face is selected. - */ - FEFacePointEvaluation( - NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, - const bool is_interior = true, - const unsigned int first_selected_component = 0); - - /** - * Reinitialize the evaluator to point to the correct precomputed mapping of - * the face in the MappingInfo object. Used in element-centric loops (ECL). - */ - void - reinit(const unsigned int cell_index, const unsigned int face_number); - - /** - * Reinitialize the evaluator to point to the correct precomputed mapping of - * the face in the MappingInfo object. Used in face-centric loops (FCL). - */ - void - reinit(const unsigned int face_index); - - /** - * This function interpolates the finite element solution, represented by - * `solution_values`, on the cell and `unit_points` passed to reinit(). - * - * @param[in] solution_values This array is supposed to contain the unknown - * values on the element read out by - * `FEEvaluation::read_dof_values(global_vector)`. - * - * @param[in] evaluation_flags Flags specifying which quantities should be - * evaluated at the points. - */ - template - void - evaluate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * This function interpolates the finite element solution, represented by - * `solution_values`, on the cell and `unit_points` passed to reinit(). - * - * @param[in] solution_values This array is supposed to contain the unknown - * values on the element as returned by `cell->get_dof_values(global_vector, - * solution_values)`. - * - * @param[in] evaluation_flags Flags specifying which quantities should be - * evaluated at the points. - */ - void - evaluate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `FEEvaluation::set_dof_values(global_vector)` or - * `FEEvaluation::distribute_local_to_global(global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - template - void - integrate(const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used to during - * `cell->set_dof_values(solution_values, global_vector)` or - * `cell->distribute_local_to_global(solution_values, global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - void - integrate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `FEEvaluation::set_dof_values(global_vector)` or - * `FEEvaluation::distribute_local_to_global(global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - template - void - test_and_sum( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used to during - * `cell->set_dof_values(solution_values, global_vector)` or - * `cell->distribute_local_to_global(solution_values, global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. Defaults to false. - * - */ - void - test_and_sum(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * Evaluate values and gradients in face for the selected face (lane) of the - * batch. Default stride into the face dofs is width of - * VectorizedArray which is the default - * vectorization over faces for FEFaceEvaluation. - */ - template - void - evaluate_in_face(const ScalarNumber *face_dof_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * Integrate values and gradients in face for the selected face (lane) of the - * batch. Default stride into the face dofs is width of - * VectorizedArray which is the default - * vectorization over faces for FEFaceEvaluation. - */ - template - void - integrate_in_face(ScalarNumber *face_dof_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * 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, Number> - normal_vector(const unsigned int point_index) const; - -private: - static constexpr std::size_t n_lanes_user_interface = - internal::VectorizedArrayTrait::width(); - static constexpr std::size_t n_lanes_internal = - internal::VectorizedArrayTrait::width(); - static constexpr std::size_t stride = - internal::VectorizedArrayTrait::stride(); - - template - void - do_evaluate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - template - void - do_integrate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); - - /** - * Actually does the evaluation templated on the chosen code path (linear or - * higher order). - */ - template - void - do_evaluate_in_face(const ScalarNumber *face_dof_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * Actually does the integration templated on the chosen code path (linear or - * higher order). - */ - template - void - do_integrate_in_face( - ScalarNumber *face_dof_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); -}; - - - template FEFacePointEvaluation:: FEFacePointEvaluation( diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index 8663c0f4f6..916e62848f 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -1577,8 +1577,10 @@ namespace NonMatching face_number.emplace_back(f_m, f_p); Assert( - cell_m->combined_face_orientation(f_m) == 1 && - cell_p->combined_face_orientation(f_p) == 1, + cell_m->combined_face_orientation(f_m) == + ReferenceCell::default_combined_face_orientation() && + cell_p->combined_face_orientation(f_p) == + ReferenceCell::default_combined_face_orientation(), ExcMessage( "Non standard face orientation is currently not implemented.")); diff --git a/tests/non_matching/mapping_info_04.cc b/tests/non_matching/mapping_info_04.cc index bee712570e..3c4297aa51 100644 --- a/tests/non_matching/mapping_info_04.cc +++ b/tests/non_matching/mapping_info_04.cc @@ -20,8 +20,6 @@ #include -#include - #include #include @@ -84,7 +82,7 @@ test_dg_fcl(const unsigned int degree, const bool curved_mesh) const unsigned int n_q_points = degree + 1; - parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + Triangulation tria; if (curved_mesh && dim > 1) GridGenerator::hyper_shell(tria, Point(), 0.5, 1, 6); diff --git a/tests/non_matching/mapping_info_05.cc b/tests/non_matching/mapping_info_05.cc index f7dc976cda..8ddcc184e9 100644 --- a/tests/non_matching/mapping_info_05.cc +++ b/tests/non_matching/mapping_info_05.cc @@ -20,8 +20,6 @@ #include -#include - #include #include @@ -80,7 +78,7 @@ test_dg_ecl(const unsigned int degree, const bool curved_mesh) const unsigned int n_q_points = degree + 1; - parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + Triangulation tria; if (curved_mesh && dim > 1) GridGenerator::hyper_shell(tria, Point(), 0.5, 1, 6); @@ -318,7 +316,7 @@ test_dg_ecl(const unsigned int degree, const bool curved_mesh) if (mask[v] == false) { for (unsigned int i = 0; - i < 2 * fe_eval_m.dofs_per_face; + i < 2 * fe_eval_m.get_dofs_projected_to_face(); ++i) fe_eval_m.get_scratch_data().begin()[i][v] = 0.; continue; -- 2.39.5