From: Peter Munch Date: Sat, 27 Nov 2021 12:01:33 +0000 (+0100) Subject: Add FEEvaluationImplHangingNodes::RunTypes::scalar X-Git-Tag: v9.4.0-rc1~786^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6b70043d56fad493f6b7bfa2b3b52fa58ca39c7f;p=dealii.git Add FEEvaluationImplHangingNodes::RunTypes::scalar --- diff --git a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h index 76cc062da0..03daf45d7f 100644 --- a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h +++ b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h @@ -38,11 +38,39 @@ namespace internal { enum class FEEvaluationImplHangingNodesRunnerTypes { + scalar, vectorized }; + /** + * Helper enum to spefify the type of vectorization for + * FEEvaluationImplHangingNodesRunnerTypes::scalar. + */ + enum class VectorizationTypes + { + /** + * Process cell by cell. + */ + index, + /** + * Process cells with the same refinement configuration together. + */ + group, + /** + * Like index but without access to individual lanes and instead use masks + * with a single entry with the value one. + */ + mask, + /** + * Assume that all lanes have the same refinement configuration. + */ + sorted + }; + + + template + struct Trait; + + template + struct Trait + { + using value_type = typename T1::value_type; + using index_type = unsigned int; + using interpolation_type = value_type; + + template + static inline const std::array, 2> & + get_interpolation_matrix(const T &fe_eval) + { + return fe_eval.get_shape_info() + .data.front() + .subface_interpolation_matrices_scalar; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + create( + const std::array mask, + const std::array + mask_new, + const unsigned int v) + { + (void)mask; + (void)mask_new; + return v; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)v; + (void)kind; + return false; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_continue(unsigned int v, + const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)v; + return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + std::array + create_mask(const std::array mask) + { + return mask; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + typename T1::value_type + get_value(const typename T1::value_type &value, const index_type &i) + { + (void)i; + return value; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + typename T1::value_type + get_value(const T1 &value, const index_type &i) + { + return value[i]; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + void + set_value(T1 & result, + const typename T1::value_type &value, + const index_type & i) + { + result[i] = value; + } + }; + + template + struct Trait + { + using value_type = T1; + using index_type = std::pair; + using interpolation_type = T1; + + template + static inline const std::array, 2> & + get_interpolation_matrix(const T &fe_eval) + { + return fe_eval.get_shape_info() + .data.front() + .subface_interpolation_matrices; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)v; + (void)kind; + return false; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_continue(unsigned int v, + const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)v; + return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + index_type + create(const std::array mask, + const std::array mask_new, + const unsigned int v) + { + (void)mask; + (void)mask_new; + T1 result = 0.0; + result[v] = 1.0; + return {result, T1(1.0) - result}; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + std::array + create_mask(const std::array mask) + { + return mask; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + T1 + get_value(const T1 &value, const index_type &) + { + return value; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + void + set_value(T1 &result, const T1 &value, const index_type &i) + { + result = result * i.second + value * i.first; + } + }; + + template + struct Trait + { + using value_type = T1; + using index_type = std::pair; + using interpolation_type = T1; + + template + static inline const std::array, 2> & + get_interpolation_matrix(const T &fe_eval) + { + return fe_eval.get_shape_info() + .data.front() + .subface_interpolation_matrices; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)v; + return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_continue(unsigned int v, + const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)v; + return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + index_type + create(const std::array mask, + const std::array mask_new, + const unsigned int v) + { + T1 result; + + for (unsigned int i = 0; i < T1::size(); ++i) + result[i] = mask_new[v] == mask[i]; + + return {result, T1(1.0) - result}; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + std::array + create_mask(const std::array mask) + { + auto new_mask = mask; + + std::sort(new_mask.begin(), new_mask.end()); + std::fill(std::unique(new_mask.begin(), new_mask.end()), + new_mask.end(), + MatrixFreeFunctions::ConstraintKinds::unconstrained); + + return new_mask; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + T1 + get_value(const T1 &value, const index_type &) + { + return value; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + void + set_value(T1 &result, const T1 &value, const index_type &i) + { + result = result * i.second + value * i.first; + } + }; + + template + struct Trait + { + using value_type = T1; + using index_type = T1; + using interpolation_type = T1; + + template + static inline const std::array, 2> & + get_interpolation_matrix(const T &fe_eval) + { + return fe_eval.get_shape_info() + .data.front() + .subface_interpolation_matrices; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)kind; + return v > 0; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + bool + do_continue(unsigned int v, + const MatrixFreeFunctions::ConstraintKinds &kind) + { + (void)kind; + + Assert(false, ExcInternalError()); + + return v > 0; // should not be called + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + T1 + create(const std::array mask, + const std::array mask_new, + const unsigned int v) + { + (void)mask; + (void)mask_new; + (void)v; + return 1.0; // return something since not used + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + std::array + create_mask(const std::array mask) + { + return mask; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + T1 + get_value(const T1 &value, const index_type &) + { + return value; + } + + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + void + set_value(T1 &result, const T1 &value, const index_type &i) + { + (void)i; + result = value; + } + }; + + + + template + class HelperBase + { + public: + inline DEAL_II_ALWAYS_INLINE + HelperBase( + const T & t, + const unsigned int & given_degree, + const bool & type_x, + const bool & type_y, + const bool & type_z, + const typename Trait::index_type &v, + const std::array< + AlignedVector< + typename Trait::interpolation_type>, + 2> & interpolation_matrices, + Number *values) + : t(t) + , given_degree(given_degree) + , type_x(type_x) + , type_y(type_y) + , type_z(type_z) + , v(v) + , interpolation_matrices(interpolation_matrices) + , values(values) + {} + + template + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + void + interpolate_3D_face( + const unsigned int dof_offset, + const unsigned int given_degree, + const typename Trait::index_type v, + const typename Trait::interpolation_type + *DEAL_II_RESTRICT weight, + Number *DEAL_II_RESTRICT values) + { + static constexpr unsigned int max_n_points_1D = 40; + + typename Trait::value_type + temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D]; + + const unsigned int points = + (fe_degree != -1 ? fe_degree : given_degree) + 1; + + AssertIndexRange(given_degree, max_n_points_1D); + + const unsigned int stride = fe_degree != -1 ? + Utilities::pow(fe_degree + 1, direction) : + Utilities::pow(given_degree + 1, direction); + + // direction side0 side1 side2 + // 0 - p^2 p + // 1 p^2 - 1 + // 2 p - 1 + const unsigned int stride2 = + ((direction == 0 && d == 1) || (direction == 1 && d == 0)) ? + (points * points) : + (((direction == 0 && d == 2) || (direction == 2 && d == 0)) ? points : + 1); + + for (unsigned int g = (skip_borders ? 1 : 0); + g < points - (skip_borders ? 1 : 0); + ++g) + { + // copy result back + for (unsigned int k = 0; k < points; ++k) + temp[k] = Trait::get_value( + values[dof_offset + k * stride + stride2 * g], v); + + // perform interpolation point by point + for (unsigned int k = 0; k < points; ++k) + { + auto sum = Trait::get_value( + weight[(transpose ? 1 : points) * k], v) * + temp[0]; + for (unsigned int h = 1; h < points; ++h) + sum += Trait::get_value( + weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h], + v) * + temp[h]; + Trait::set_value( + values[dof_offset + k * stride + stride2 * g], sum, v); + } + } + } + + template + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + void + interpolate_3D_edge( + const unsigned int p, + const unsigned int given_degree, + const typename Trait::index_type v, + const typename Trait::interpolation_type + *DEAL_II_RESTRICT weight, + Number *DEAL_II_RESTRICT values) + { + static constexpr unsigned int max_n_points_1D = 40; + + typename Trait::value_type + temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D]; + + const unsigned int points = + (fe_degree != -1 ? fe_degree : given_degree) + 1; + + AssertIndexRange(given_degree, max_n_points_1D); + + const unsigned int stride = fe_degree != -1 ? + Utilities::pow(fe_degree + 1, direction) : + Utilities::pow(given_degree + 1, direction); + + // copy result back + for (unsigned int k = 0; k < points; ++k) + temp[k] = + Trait::get_value(values[p + k * stride], + v); + + // perform interpolation point by point + for (unsigned int k = 0; k < points; ++k) + { + auto sum = Trait::get_value( + weight[(transpose ? 1 : points) * k], v) * + temp[0]; + for (unsigned int h = 1; h < points; ++h) + sum += Trait::get_value( + weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h], + v) * + temp[h]; + Trait::set_value(values[p + k * stride], + sum, + v); + } + } + + template + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + void + process_edge() const + { + if (do_x) + interpolate_3D_edge<0>(t.line(0, type_y, type_z), + given_degree, + v, + interpolation_matrices[!type_x].data(), + values); + + if (do_y) + interpolate_3D_edge<1>(t.line(1, type_x, type_z), + given_degree, + v, + interpolation_matrices[!type_y].data(), + values); + + if (do_z) + interpolate_3D_edge<2>(t.line(2, type_x, type_y), + given_degree, + v, + interpolation_matrices[!type_z].data(), + values); + } + + template + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + void + process_faces_fast() const + { + static_assert((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) || + (!do_x && !do_y && do_z), + "Only one face can be chosen."); + + static const unsigned int direction = do_x ? 0 : (do_y ? 1 : 2); + const bool type = do_x ? type_x : (do_y ? type_y : type_z); + + if (!do_x) + interpolate_3D_face<0, direction, false>( + t.face(direction, type), + given_degree, + v, + interpolation_matrices[!type_x].data(), + values); + + if (!do_y) + interpolate_3D_face<1, direction, false>( + t.face(direction, type), + given_degree, + v, + interpolation_matrices[!type_y].data(), + values); + + if (!do_z) + interpolate_3D_face<2, direction, false>( + t.face(direction, type), + given_degree, + v, + interpolation_matrices[!type_z].data(), + values); + } + + template + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + void + process_faces() const + { + static_assert(((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) || + (!do_x && !do_y && do_z)) == false, + "Only one face can be chosen."); + + // direction 0 + { + const auto inpterolation_matrix = + interpolation_matrices[!type_x].data(); + + // faces + if (do_y && given_degree > 1) + interpolate_3D_face<0, 1, true>( + t.face(1, type_y), given_degree, v, inpterolation_matrix, values); + + if (do_z && given_degree > 1) + interpolate_3D_face<0, 2, true>( + t.face(2, type_z), given_degree, v, inpterolation_matrix, values); + + // direction 0 -> edges + interpolate_3D_edge<0>((do_x && do_y && !do_z) ? + (t.lines_plane(0, type_x, type_y, 0)) : + ((do_x && !do_y && do_z) ? + (t.lines_plane(1, type_x, type_z, 0)) : + (t.lines(0, type_y, type_z, 0))), + given_degree, + v, + inpterolation_matrix, + values); + + + interpolate_3D_edge<0>((do_x && do_y && !do_z) ? + (t.lines_plane(0, type_x, type_y, 1)) : + ((do_x && !do_y && do_z) ? + (t.lines_plane(1, type_x, type_z, 1)) : + (t.lines(0, type_y, type_z, 1))), + given_degree, + v, + inpterolation_matrix, + values); + + if (do_y && do_z) + interpolate_3D_edge<0>(t.lines(0, type_y, type_z, 2), + given_degree, + v, + inpterolation_matrix, + values); + } + + // direction 1 + { + const auto inpterolation_matrix = + interpolation_matrices[!type_y].data(); + + // faces + if (do_x && given_degree > 1) + interpolate_3D_face<1, 0, true>( + t.face(0, type_x), given_degree, v, inpterolation_matrix, values); + + if (do_z && given_degree > 1) + interpolate_3D_face<1, 2, true>( + t.face(2, type_z), given_degree, v, inpterolation_matrix, values); + + // lines + interpolate_3D_edge<1>((do_x && do_y && !do_z) ? + (t.lines_plane(0, type_x, type_y, 2)) : + ((!do_x && do_y && do_z) ? + (t.lines_plane(2, type_y, type_z, 0)) : + (t.lines(1, type_x, type_z, 0))), + given_degree, + v, + inpterolation_matrix, + values); + + interpolate_3D_edge<1>((do_x && do_y && !do_z) ? + (t.lines_plane(0, type_x, type_y, 3)) : + ((!do_x && do_y && do_z) ? + (t.lines_plane(2, type_y, type_z, 1)) : + (t.lines(1, type_x, type_z, 1))), + given_degree, + v, + inpterolation_matrix, + values); + + if (do_x && do_z) + interpolate_3D_edge<1>(t.lines(1, type_x, type_z, 2), + given_degree, + v, + inpterolation_matrix, + values); + } + + // direction 2 -> faces + { + const auto inpterolation_matrix = + interpolation_matrices[!type_z].data(); + + if (do_x && given_degree > 1) + interpolate_3D_face<2, 0, true>( + t.face(0, type_x), given_degree, v, inpterolation_matrix, values); + + if (do_y && given_degree > 1) + interpolate_3D_face<2, 1, true>( + t.face(1, type_y), given_degree, v, inpterolation_matrix, values); + + // direction 2 -> edges + interpolate_3D_edge<2>((do_x && !do_y && do_z) ? + (t.lines_plane(1, type_x, type_z, 2)) : + ((!do_x && do_y && do_z) ? + (t.lines_plane(2, type_y, type_z, 2)) : + (t.lines(2, type_x, type_y, 0))), + given_degree, + v, + inpterolation_matrix, + values); + + interpolate_3D_edge<2>((do_x && !do_y && do_z) ? + (t.lines_plane(1, type_x, type_z, 3)) : + ((!do_x && do_y && do_z) ? + (t.lines_plane(2, type_y, type_z, 3)) : + (t.lines(2, type_x, type_y, 1))), + given_degree, + v, + inpterolation_matrix, + values); + + if (do_x && do_y) + interpolate_3D_edge<2>(t.lines(2, type_x, type_y, 2), + given_degree, + v, + inpterolation_matrix, + values); + } + } + + private: + const T & t; + const unsigned int & given_degree; + const bool & type_x; + const bool & type_y; + const bool & type_z; + const typename Trait::index_type &v; + const std::array< + AlignedVector< + typename Trait::interpolation_type>, + 2> & interpolation_matrices; + Number *values; + }; + + /** + * Helper enum to specify which Helper implementation should be used. + */ + enum class HelperType + { + /** + * Compute the start indices of faces and edges based on the template + * argument fe_degree. + */ + constant, + /** + * Compute the start indices of faces and edges based on the fe_degree + * passed to the constructor (to be used if the template argument is -1). + */ + dynamic + }; + + template + class Helper; + + template + class Helper : public HelperBase, + Number, + VectorizationType, + fe_degree, + transpose> + { + public: +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + Helper(const unsigned int &given_degree, + const bool & type_x, + const bool & type_y, + const bool & type_z, + const typename Trait::index_type &v, + const std::array< + AlignedVector< + typename Trait::interpolation_type>, + 2> & interpolation_matrices, + Number *values) + : HelperBase, + Number, + VectorizationType, + fe_degree, + transpose>(*this, + given_degree, + type_x, + type_y, + type_z, + v, + interpolation_matrices, + values) + , points(given_degree + 1) + { + static_assert(fe_degree == -1, "Only working for fe_degree = -1."); + } + + const unsigned int points; + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + line(unsigned int i, unsigned int j, unsigned int k) const + { + return line_array[i][j][k]; + } + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + face(unsigned int i, unsigned int j) const + { + return face_array[i][j]; + } + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + lines_plane(unsigned int i, + unsigned int j, + unsigned int k, + unsigned int l) const + { + return lines_plane_array[i][j][k][l]; + } + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + lines(unsigned int i, + unsigned int j, + unsigned int k, + unsigned int l) const + { + return lines_array[i][j][k][l]; + } + + private: + const dealii::ndarray line_array = { + {{{{{points * points * points - points, points *points - points}}, + {{points * points * points - points * points, 0}}}}, + {{{{points * points * points - points * points + points - 1, + points - 1}}, + {{points * points * points - points * points, 0}}}}, + {{{{points * points - 1, points - 1}}, + {{points * points - points, 0}}}}}}; + + const dealii::ndarray face_array = { + {{{points - 1, 0}}, + {{points * points - points, 0}}, + {{points * points * points - points * points, 0}}}}; + + const dealii::ndarray lines_plane_array = { + {{{{{{{points * points - points, + points *points *points - points, + points - 1, + points *points *points - points *points + points - 1}}, + {{0, + points *points *points - points *points, + points - 1, + points *points *points - points *points + points - 1}}}}, + {{{{points * points - points, + points *points *points - points, + 0, + points *points *points - points *points}}, + {{0, + points *points *points - points *points, + 0, + points *points *points - points *points}}}}}}, + {{{{{{points * points * points - points * points, + points *points *points - points, + points - 1, + points *points - 1}}, + {{0, points *points - points, points - 1, points *points - 1}}}}, + {{{{points * points * points - points * points, + points *points *points - points, + 0, + points *points - points}}, + {{0, points *points - points, 0, points *points - points}}}}}}, + {{{{{{points * points * points - points * points, + points *points *points - points *points + points - 1, + points * points - points, + points * points - 1}}, + {{0, points - 1, points *points - points, points *points - 1}}}}, + {{{{points * points * points - points * points, + points *points *points - points *points + points - 1, + 0, + points - 1}}, + {{0, points - 1, 0, points - 1}}}}}}}}; + + const dealii::ndarray lines_array = { + {{{{{{{points * points - points, + points *points *points - points *points, + points *points *points - points}}, + {{0, points *points - points, points *points *points - points}}}}, + {{{{0, + points *points *points - points *points, + points *points *points - points}}, + {{0, + points *points - points, + points *points *points - points *points}}}}}}, + {{{{{{points - 1, + points *points *points - points *points, + points *points *points - points *points + points - 1}}, + {{0, + points - 1, + points *points *points - points *points + points - 1}}}}, + {{{{0, + points *points *points - points *points, + points *points *points - points *points + points - 1}}, + {{0, points - 1, points *points *points - points *points}}}}}}, + {{{{{{points - 1, points *points - points, points *points - 1}}, + {{0, points - 1, points *points - 1}}}}, + {{{{0, points *points - points, points *points - 1}}, + {{0, points - 1, points *points - points}}}}}}}}; + }; + + template + class Helper : public HelperBase, + Number, + VectorizationType, + fe_degree, + transpose> + { + public: +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + Helper(const unsigned int &given_degree, + const bool & type_x, + const bool & type_y, + const bool & type_z, + const typename Trait::index_type &v, + const std::array< + AlignedVector< + typename Trait::interpolation_type>, + 2> & interpolation_matrices, + Number *values) + : HelperBase, + Number, + VectorizationType, + fe_degree, + transpose>(*this, + given_degree, + type_x, + type_y, + type_z, + v, + interpolation_matrices, + values) + { + static_assert(fe_degree != -1, "Only working for fe_degree != -1."); + } + + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + line(unsigned int i, unsigned int j, unsigned int k) const + { + static constexpr unsigned int points = fe_degree + 1; + + static constexpr dealii::ndarray line_array = { + {{{{{points * points * points - points, points * points - points}}, + {{points * points * points - points * points, 0}}}}, + {{{{points * points * points - points * points + points - 1, + points - 1}}, + {{points * points * points - points * points, 0}}}}, + {{{{points * points - 1, points - 1}}, + {{points * points - points, 0}}}}}}; + + return line_array[i][j][k]; + } + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + face(unsigned int i, unsigned int j) const + { + static constexpr unsigned int points = fe_degree + 1; + + static constexpr dealii::ndarray face_array = { + {{{points - 1, 0}}, + {{points * points - points, 0}}, + {{points * points * points - points * points, 0}}}}; + + return face_array[i][j]; + } + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + lines_plane(unsigned int i, + unsigned int j, + unsigned int k, + unsigned int l) const + { + static constexpr unsigned int points = fe_degree + 1; + + static constexpr dealii::ndarray + lines_plane_array = { + {{{{{{{points * points - points, + points * points * points - points, + points - 1, + points * points * points - points * points + points - 1}}, + {{0, + points * points * points - points * points, + points - 1, + points * points * points - points * points + points - 1}}}}, + {{{{points * points - points, + points * points * points - points, + 0, + points * points * points - points * points}}, + {{0, + points * points * points - points * points, + 0, + points * points * points - points * points}}}}}}, + {{{{{{points * points * points - points * points, + points * points * points - points, + points - 1, + points * points - 1}}, + {{0, + points * points - points, + points - 1, + points * points - 1}}}}, + {{{{points * points * points - points * points, + points * points * points - points, + 0, + points * points - points}}, + {{0, points * points - points, 0, points * points - points}}}}}}, + {{{{{{points * points * points - points * points, + points * points * points - points * points + points - 1, + points * points - points, + points * points - 1}}, + {{0, + points - 1, + points * points - points, + points * points - 1}}}}, + {{{{points * points * points - points * points, + points * points * points - points * points + points - 1, + 0, + points - 1}}, + {{0, points - 1, 0, points - 1}}}}}}}}; + + return lines_plane_array[i][j][k][l]; + } + +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + unsigned int + lines(unsigned int i, + unsigned int j, + unsigned int k, + unsigned int l) const + { + static constexpr unsigned int points = fe_degree + 1; + + static constexpr dealii::ndarray lines_array = { + {{{{{{{points * points - points, + points * points * points - points * points, + points * points * points - points}}, + {{0, + points * points - points, + points * points * points - points}}}}, + {{{{0, + points * points * points - points * points, + points * points * points - points}}, + {{0, + points * points - points, + points * points * points - points * points}}}}}}, + {{{{{{points - 1, + points * points * points - points * points, + points * points * points - points * points + points - 1}}, + {{0, + points - 1, + points * points * points - points * points + points - 1}}}}, + {{{{0, + points * points * points - points * points, + points * points * points - points * points + points - 1}}, + {{0, points - 1, points * points * points - points * points}}}}}}, + {{{{{{points - 1, points * points - points, points * points - 1}}, + {{0, points - 1, points * points - 1}}}}, + {{{{0, points * points - points, points * points - 1}}, + {{0, points - 1, points * points - points}}}}}}}}; + + return lines_array[i][j][k][l]; + } + }; + + + template + class FEEvaluationImplHangingNodesRunner< + FEEvaluationImplHangingNodesRunnerTypes::scalar, + dim, + fe_degree, + Number, + is_face> + { + public: + static const VectorizationTypes VectorizationType = + VectorizationTypes::index; + + private: + template + static +#ifndef DEBUG + inline DEAL_II_ALWAYS_INLINE +#endif + void + interpolate_2D( + const unsigned int given_degree, + const typename Trait::index_type v, + const typename Trait::interpolation_type + *DEAL_II_RESTRICT weight, + Number *DEAL_II_RESTRICT values) + { + static constexpr unsigned int max_n_points_1D = 40; + + typename Trait::value_type + temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D]; + + const unsigned int points = + (fe_degree != -1 ? fe_degree : given_degree) + 1; + + AssertIndexRange(given_degree, max_n_points_1D); + + const unsigned int d = side / 2; // direction + const unsigned int s = side % 2; // left or right surface + + const unsigned int offset = dealii::Utilities::pow(points, d + 1); + const unsigned int stride = + (s == 0 ? 0 : (points - 1)) * dealii::Utilities::pow(points, d); + + const unsigned int r1 = dealii::Utilities::pow(points, dim - d - 1); + const unsigned int r2 = dealii::Utilities::pow(points, d); + + // copy result back + for (unsigned int i = 0, k = 0; i < r1; ++i) + for (unsigned int j = 0; j < r2; ++j, ++k) + temp[k] = Trait::get_value( + values[i * offset + stride + j], v); + + // perform interpolation point by point (note: r1 * r2 == + // points^(dim-1)) + for (unsigned int i = 0, k = 0; i < r1; ++i) + for (unsigned int j = 0; j < r2; ++j, ++k) + { + typename Trait::value_type sum = 0.0; + for (unsigned int h = 0; h < points; ++h) + sum += Trait::get_value( + weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h], + v) * + temp[h]; + Trait::set_value( + values[i * offset + stride + j], sum, v); + } + } + + public: + template + static void + run_internal(const unsigned int n_desired_components, + const FEEvaluationBaseData &fe_eval, + const std::array & constraint_mask, + Number * values) + { + const unsigned int given_degree = + fe_degree != -1 ? fe_degree : + fe_eval.get_shape_info().data.front().fe_degree; + + const auto &interpolation_matrices = + Trait::get_interpolation_matrix(fe_eval); + + const auto constraint_mask_sorted = + Trait::create_mask(constraint_mask); + + for (unsigned int c = 0; c < n_desired_components; ++c) + { + for (unsigned int v = 0; v < Number::size(); ++v) + { + const auto mask = constraint_mask_sorted[v]; + + if (Trait::do_break(v, mask)) + break; + + if (Trait::do_continue(v, mask)) + continue; + + const auto vv = + Trait::create(constraint_mask, + constraint_mask_sorted, + v); + + if (dim == 2) // 2D: only faces + { + const auto is_set = [](const auto a, const auto b) -> bool { + return (a & b) == b; + }; + + // direction 0: + if ((mask & MatrixFreeFunctions::ConstraintKinds::face_y) != + MatrixFreeFunctions::ConstraintKinds::unconstrained) + { + const bool is_subface_0 = + (mask & + MatrixFreeFunctions::ConstraintKinds::subcell_x) == + MatrixFreeFunctions::ConstraintKinds::unconstrained; + + const auto *weights = + interpolation_matrices[is_subface_0].data(); + + if (is_set( + mask, + MatrixFreeFunctions::ConstraintKinds::subcell_y)) + interpolate_2D<2, transpose>(given_degree, + vv, + weights, + values); // face 2 + else + interpolate_2D<3, transpose>(given_degree, + vv, + weights, + values); // face 3 + } + + // direction 1: + if ((mask & MatrixFreeFunctions::ConstraintKinds::face_x) != + MatrixFreeFunctions::ConstraintKinds::unconstrained) + { + const bool is_subface_0 = + (mask & + MatrixFreeFunctions::ConstraintKinds::subcell_y) == + MatrixFreeFunctions::ConstraintKinds::unconstrained; + + const auto *weights = + interpolation_matrices[is_subface_0].data(); + + if (is_set( + mask, + MatrixFreeFunctions::ConstraintKinds::subcell_x)) + interpolate_2D<0, transpose>(given_degree, + vv, + weights, + values); // face 0 + else + interpolate_2D<1, transpose>(given_degree, + vv, + weights, + values); // face 1 + } + } + else if (dim == 3) // 3D faces and edges + { + const auto m = static_cast(mask); + + const bool type_x = (m >> 0) & 1; + const bool type_y = (m >> 1) & 1; + const bool type_z = (m >> 2) & 1; + const auto faces = (m >> 3) & 7; + const auto edges = (m >> 6); + + Helper + helper(given_degree, + type_x, + type_y, + type_z, + vv, + interpolation_matrices, + values); + + if (faces > 0) + switch (faces) + { + case 0: + break; + case 1: + helper + .template process_faces_fast(); + break; + case 2: + helper + .template process_faces_fast(); + break; + case 3: + helper.template process_faces(); + break; + case 4: + helper + .template process_faces_fast(); + break; + case 5: + helper.template process_faces(); + break; + case 6: + helper.template process_faces(); + break; + case 7: + helper.template process_faces(); + break; + } + + if (edges > 0) + switch (edges) + { + case 0: + break; + case 1: + helper.template process_edge(); + break; + case 2: + helper.template process_edge(); + break; + case 3: + helper.template process_edge(); + break; + case 4: + helper.template process_edge(); + break; + case 5: + helper.template process_edge(); + break; + case 6: + helper.template process_edge(); + break; + case 7: + helper.template process_edge(); + break; + } + } + else + { + Assert(false, ExcNotImplemented()); + } + } + + values += fe_eval.get_shape_info().dofs_per_component_on_cell; + } + } + }; + template @@ -449,12 +1891,12 @@ namespace internal & c_mask, Number *values) { - using RunnerType = FEEvaluationImplHangingNodesRunner< - FEEvaluationImplHangingNodesRunnerTypes::vectorized, - dim, - fe_degree, - Number, - is_face>; + using RunnerType = + FEEvaluationImplHangingNodesRunner(), + dim, + fe_degree, + Number, + is_face>; if (transpose) RunnerType::template run_internal(n_desired_components, @@ -469,6 +1911,15 @@ namespace internal return false; } + + template + static constexpr FEEvaluationImplHangingNodesRunnerTypes + used_runner_type() + { + return ((Number::size() > 2) && (fe_degree == -1 || fe_degree > 2)) ? + FEEvaluationImplHangingNodesRunnerTypes::vectorized : + FEEvaluationImplHangingNodesRunnerTypes::scalar; + } };