const bool is_interior;
};
-// ----------------------- template and inline function ----------------------
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
- FEPointEvaluationBase(const Mapping<dim> &mapping,
- const FiniteElement<dim> &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<NonMatching::MappingInfo<dim, spacedim, Number>>(
- 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 <int n_components_,
+ int dim,
+ int spacedim = dim,
+ typename Number = double>
+class FEPointEvaluation
+ : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
{
- 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<Number>::value_type;
+ using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
+ using ETT = typename internal::FEPointEvaluation::
+ EvaluatorTypeTraits<dim, n_components, Number>;
+ 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 <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
- 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<dim> &mapping,
+ const FiniteElement<dim> &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<dim, spacedim, Number> &mapping_info,
const FiniteElement<dim> &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<dim, spacedim>::cell_iterator &cell,
+ const ArrayView<const Point<dim>> &unit_points);
+ /**
+ * Reinitialize the evaluator to point to the correct precomputed mapping of
+ * the single cell in the MappingInfo object.
+ */
+ void
+ reinit();
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
- FEPointEvaluationBase(
- FEPointEvaluationBase<n_components_, dim, spacedim, Number> &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<NonMatching::MappingInfo<dim, spacedim, Number>>(
- *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 <std::size_t stride_view>
+ void
+ evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
- FEPointEvaluationBase(
- FEPointEvaluationBase<n_components_, dim, spacedim, Number>
- &&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<const ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
- ~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 <std::size_t stride_view>
+ void
+ integrate(const StridedArrayView<ScalarNumber, stride_view> &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<ScalarNumber> &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 <std::size_t stride_view>
+ void
+ test_and_sum(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
-template <int n_components_, int dim, int spacedim, typename Number>
-void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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<ScalarNumber> &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<Number>::width();
+ static constexpr std::size_t n_lanes_internal =
+ internal::VectorizedArrayTrait<VectorizedArrayType>::width();
+ static constexpr std::size_t stride =
+ internal::VectorizedArrayTrait<Number>::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 <bool is_linear, std::size_t stride_view>
+ void
+ prepare_evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &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 <bool is_linear, std::size_t stride_view>
+ void
+ compute_evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &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 <bool is_linear, std::size_t stride_view>
+ void
+ evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &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 <std::size_t stride_view>
+ void
+ evaluate_slow(
+ const StridedArrayView<const ScalarNumber, stride_view> &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 <bool is_linear>
+ 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 <bool is_linear, std::size_t stride_view>
+ void
+ finish_integrate_fast(
+ const StridedArrayView<ScalarNumber, stride_view> &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 <bool do_JxW, bool is_linear, std::size_t stride_view>
+ void
+ integrate_fast(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
- fast_path = false;
- }
-}
+ /**
+ * Slow path of the integrate function using FEValues.
+ */
+ template <bool do_JxW, std::size_t stride_view>
+ void
+ integrate_slow(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+ /**
+ * Implementation of the integrate/test_and_sum function.
+ */
+ template <bool do_JxW, std::size_t stride_view>
+ void
+ do_integrate(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+};
-template <int n_components_, int dim, int spacedim, typename Number>
-template <bool is_face, bool is_linear>
-inline void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::do_reinit()
-{
- const unsigned int geometry_index =
- mapping_info->template compute_geometry_index_offset<is_face>(
- 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 <int n_components_,
+ int dim,
+ int spacedim = dim,
+ typename Number = double>
+class FEFacePointEvaluation
+ : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
+{
+public:
+ static constexpr unsigned int dimension = dim;
+ static constexpr unsigned int n_components = n_components_;
- const_cast<unsigned int &>(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<unsigned int &>(n_q_batches) =
- (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
+ using ScalarNumber =
+ typename internal::VectorizedArrayTrait<Number>::value_type;
+ using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
+ using ETT = typename internal::FEPointEvaluation::
+ EvaluatorTypeTraits<dim, n_components, Number>;
+ 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<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> &fe,
+ const bool is_interior = true,
+ const unsigned int first_selected_component = 0);
- const_cast<unsigned int &>(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 <std::size_t stride_view>
+ void
+ evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &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<const ScalarNumber> &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 <std::size_t stride_view>
+ void
+ integrate(const StridedArrayView<ScalarNumber, stride_view> &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<ScalarNumber> &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 <int n_components_, int dim, int spacedim, typename Number>
-inline const typename FEPointEvaluationBase<n_components_,
- dim,
- spacedim,
- Number>::value_type &
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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 <std::size_t stride_view>
+ void
+ test_and_sum(
+ const StridedArrayView<ScalarNumber, stride_view> &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<ScalarNumber> &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<selected_floating_point_type> which is the default
+ * vectorization over faces for FEFaceEvaluation.
+ */
+ template <int stride_face_dof = VectorizedArrayType::size()>
+ void
+ evaluate_in_face(const ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
-template <int n_components_, int dim, int spacedim, typename Number>
-inline const typename FEPointEvaluationBase<n_components_,
- dim,
- spacedim,
- Number>::gradient_type &
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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<selected_floating_point_type> which is the default
+ * vectorization over faces for FEFaceEvaluation.
+ */
+ template <int stride_face_dof = VectorizedArrayType::size()>
+ 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<Number>::width();
+ static constexpr std::size_t n_lanes_internal =
+ internal::VectorizedArrayTrait<VectorizedArrayType>::width();
+ static constexpr std::size_t stride =
+ internal::VectorizedArrayTrait<Number>::stride();
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_value(
- const value_type &value,
- const unsigned int point_index)
-{
- AssertIndexRange(point_index, n_q_points);
- values[point_index] = value;
-}
+ template <bool is_linear, std::size_t stride_view>
+ void
+ do_evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+ template <bool do_JxW, bool is_linear, std::size_t stride_view>
+ void
+ do_integrate(
+ const StridedArrayView<ScalarNumber, stride_view> &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 <bool is_linear, int stride_face_dof>
+ void
+ do_evaluate_in_face(const ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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 <bool do_JxW, bool is_linear, int stride_face_dof>
+ void
+ do_integrate_in_face(
+ ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+};
+// ----------------------- template and inline function ----------------------
template <int n_components_, int dim, int spacedim, typename Number>
-inline DerivativeForm<1, dim, spacedim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::jacobian(
- const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(const Mapping<dim> &mapping,
+ const FiniteElement<dim> &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<NonMatching::MappingInfo<dim, spacedim, Number>>(
+ 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 <int n_components_, int dim, int spacedim, typename Number>
-inline DerivativeForm<1, spacedim, dim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::inverse_jacobian(
- const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> &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 <int n_components_, int dim, int spacedim, typename Number>
-inline Number
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::JxW(
- const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(
+ FEPointEvaluationBase<n_components_, dim, spacedim, Number> &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<NonMatching::MappingInfo<dim, spacedim, Number>>(
+ *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 <int n_components_, int dim, int spacedim, typename Number>
-inline Point<spacedim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::real_point(
- const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(
+ FEPointEvaluationBase<n_components_, dim, spacedim, Number>
+ &&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 <int n_components_, int dim, int spacedim, typename Number>
-inline Point<dim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::unit_point(
- const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ ~FEPointEvaluationBase()
{
- AssertIndexRange(point_index, n_q_points);
- Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
- Point<dim, Number> unit_point;
- for (unsigned int d = 0; d < dim; ++d)
- unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
- unit_point_ptr[point_index / stride][d], point_index % stride);
- return unit_point;
+ connection_is_reinitialized.disconnect();
}
template <int n_components_, int dim, int spacedim, typename Number>
-inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
- quadrature_point_indices() const
+void
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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 <int n_components_,
- int dim,
- int spacedim = dim,
- typename Number = double>
-class FEPointEvaluation
- : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
-{
-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<Number>::value_type;
- using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
- Number>::vectorized_value_type;
- using ETT = typename internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>;
- 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<dim> &mapping,
- const FiniteElement<dim> &fe,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component = 0)
- : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
- 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<dim, spacedim, Number> &mapping_info,
- const FiniteElement<dim> &fe,
- const unsigned int first_selected_component = 0)
- : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
- 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<dim, spacedim>::cell_iterator &cell,
- const ArrayView<const Point<dim>> &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 <std::size_t stride_view>
- void
- evaluate(
- const StridedArrayView<const ScalarNumber, stride_view> &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<const ScalarNumber> &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 <std::size_t stride_view>
- void
- integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values = false);
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face, bool is_linear>
+inline void
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::do_reinit()
+{
+ const unsigned int geometry_index =
+ mapping_info->template compute_geometry_index_offset<is_face>(
+ 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<ScalarNumber> &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 <std::size_t stride_view>
- void
- test_and_sum(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values = false);
+ const_cast<unsigned int &>(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<ScalarNumber> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values = false);
+ // round up n_q_points_scalar / n_lanes_internal
+ const_cast<unsigned int &>(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<Number>::width();
- static constexpr std::size_t n_lanes_internal =
- internal::VectorizedArrayTrait<VectorizedArrayType>::width();
- static constexpr std::size_t stride =
- internal::VectorizedArrayTrait<Number>::stride();
+ const_cast<unsigned int &>(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 <bool is_linear, std::size_t stride_view>
- void
- prepare_evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &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 <bool is_linear, std::size_t stride_view>
- void
- compute_evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &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 <bool is_linear, std::size_t stride_view>
- void
- evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &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 <std::size_t stride_view>
- void
- evaluate_slow(
- const StridedArrayView<const ScalarNumber, stride_view> &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 <bool is_linear>
- 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 <bool is_linear, std::size_t stride_view>
- void
- finish_integrate_fast(
- const StridedArrayView<ScalarNumber, stride_view> &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 <bool do_JxW, bool is_linear, std::size_t stride_view>
- void
- integrate_fast(
- const StridedArrayView<ScalarNumber, stride_view> &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 <bool do_JxW, std::size_t stride_view>
- void
- integrate_slow(
- const StridedArrayView<ScalarNumber, stride_view> &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 <bool do_JxW, std::size_t stride_view>
- void
- do_integrate(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values);
-};
+ is_reinitialized = true;
+}
template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit()
+inline const typename FEPointEvaluationBase<n_components_,
+ dim,
+ spacedim,
+ Number>::value_type &
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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<false, true>();
- else
- this->template do_reinit<false, false>();
+ AssertIndexRange(point_index, values.size());
+ return values[point_index];
}
template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ArrayView<const Point<dim>> &unit_points)
+inline const typename FEPointEvaluationBase<n_components_,
+ dim,
+ spacedim,
+ Number>::gradient_type &
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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<FEValues<dim, spacedim>>(
- *this->mapping,
- *this->fe,
- Quadrature<dim>(
- std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
- this->update_flags);
- this->fe_values->reinit(cell);
- }
- if (this->use_linear_path)
- this->template do_reinit<false, true>();
- else
- this->template do_reinit<false, false>();
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_value(
+ const value_type &value,
+ const unsigned int point_index)
+{
+ AssertIndexRange(point_index, n_q_points);
+ values[point_index] = value;
}
template <int n_components_, int dim, int spacedim, typename Number>
inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
- const unsigned int cell_index)
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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<false, true>();
- else
- this->template do_reinit<false, false>();
-
- if (!this->fast_path)
- {
- std::vector<Point<dim>> 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<FEValues<dim, spacedim>>(
- *this->mapping,
- *this->fe,
- Quadrature<dim>(
- std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
- this->update_flags);
- this->fe_values->reinit(
- this->mapping_info->get_cell_iterator(this->current_cell_index));
- }
+template <int n_components_, int dim, int spacedim, typename Number>
+inline DerivativeForm<1, dim, spacedim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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 <int n_components_, int dim, int spacedim, typename Number>
-template <std::size_t stride_view>
-void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags)
+inline DerivativeForm<1, spacedim, dim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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 <int n_components_, int dim, int spacedim, typename Number>
+inline Number
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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 <int n_components_, int dim, int spacedim, typename Number>
+inline Point<spacedim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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 <int n_components_, int dim, int spacedim, typename Number>
+inline Point<dim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::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<dim, Number> unit_point;
+ for (unsigned int d = 0; d < dim; ++d)
+ unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+ unit_point_ptr[point_index / stride][d], point_index % stride);
+ return unit_point;
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ quadrature_point_indices() const
+{
+ return {0U, n_q_points};
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> &fe,
+ const unsigned int first_selected_component)
+ : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
+ mapping_info,
+ fe,
+ first_selected_component)
+{}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
+ const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component)
+ : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
+ mapping,
+ fe,
+ update_flags,
+ first_selected_component)
+{}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::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<false, true>();
+ else
+ this->template do_reinit<false, false>();
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const ArrayView<const Point<dim>> &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<FEValues<dim, spacedim>>(
+ *this->mapping,
+ *this->fe,
+ Quadrature<dim>(
+ std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+ this->update_flags);
+ this->fe_values->reinit(cell);
+ }
+
+ if (this->use_linear_path)
+ this->template do_reinit<false, true>();
+ else
+ this->template do_reinit<false, false>();
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::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<false, true>();
+ else
+ this->template do_reinit<false, false>();
+
+ if (!this->fast_path)
+ {
+ std::vector<Point<dim>> 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<FEValues<dim, spacedim>>(
+ *this->mapping,
+ *this->fe,
+ Quadrature<dim>(
+ std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+ this->update_flags);
+
+ this->fe_values->reinit(
+ this->mapping_info->get_cell_iterator(this->current_cell_index));
+ }
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &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<true>(solution_values, evaluation_flags);
-/**
- * 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 <int n_components_,
- int dim,
- int spacedim = dim,
- typename Number = double>
-class FEFacePointEvaluation
- : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
-{
-public:
- static constexpr unsigned int dimension = dim;
- static constexpr unsigned int n_components = n_components_;
-
- using number_type = Number;
-
- using ScalarNumber =
- typename internal::VectorizedArrayTrait<Number>::value_type;
- using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
- Number>::vectorized_value_type;
- using ETT = typename internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>;
- 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<dim, spacedim, Number> &mapping_info,
- const FiniteElement<dim> &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 <std::size_t stride_view>
- void
- evaluate(
- const StridedArrayView<const ScalarNumber, stride_view> &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<const ScalarNumber> &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 <std::size_t stride_view>
- void
- integrate(const StridedArrayView<ScalarNumber, stride_view> &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<ScalarNumber> &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 <std::size_t stride_view>
- void
- test_and_sum(
- const StridedArrayView<ScalarNumber, stride_view> &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<ScalarNumber> &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<selected_floating_point_type> which is the default
- * vectorization over faces for FEFaceEvaluation.
- */
- template <int stride_face_dof = VectorizedArrayType::size()>
- 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<selected_floating_point_type> which is the default
- * vectorization over faces for FEFaceEvaluation.
- */
- template <int stride_face_dof = VectorizedArrayType::size()>
- 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<Number>::width();
- static constexpr std::size_t n_lanes_internal =
- internal::VectorizedArrayTrait<VectorizedArrayType>::width();
- static constexpr std::size_t stride =
- internal::VectorizedArrayTrait<Number>::stride();
-
- template <bool is_linear, std::size_t stride_view>
- void
- do_evaluate(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags);
-
- template <bool do_JxW, bool is_linear, std::size_t stride_view>
- void
- do_integrate(
- const StridedArrayView<ScalarNumber, stride_view> &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 <bool is_linear, int stride_face_dof>
- 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 <bool do_JxW, bool is_linear, int stride_face_dof>
- void
- do_integrate_in_face(
- ScalarNumber *face_dof_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values);
-};
-
-
-
template <int n_components_, int dim, int spacedim, typename Number>
FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
FEFacePointEvaluation(
projects / dealii.git / commitdiff