From: Peter Munch <peterrmuench@gmail.com> Date: Tue, 28 Sep 2021 10:51:56 +0000 (+0200) Subject: Split up evaluation_template_factory X-Git-Tag: v9.4.0-rc1~974^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12778%2Fhead;p=dealii.git Split up evaluation_template_factory --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 2fa379d2b9..215a3ecd48 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -4906,587 +4906,6 @@ namespace internal } }; - - - template <int dim, typename Number, bool is_face> - struct FEEvaluationImplHangingNodes - { - template <int fe_degree, int n_q_points_1d> - static bool - run(const unsigned int n_desired_components, - const FEEvaluationBaseData<dim, - typename Number::value_type, - is_face, - Number> &fe_eval, - const bool transpose, - const std::array<MatrixFreeFunctions::ConstraintKinds, Number::size()> - & c_mask, - Number *values) - { - if (transpose) - run_internal<fe_degree, true>(n_desired_components, - fe_eval, - c_mask, - values); - else - run_internal<fe_degree, false>(n_desired_components, - fe_eval, - c_mask, - values); - - return false; - } - - private: - template <int fe_degree, unsigned int side, bool transpose> - static void - interpolate_2D(const unsigned int given_degree, - const unsigned int v, - const Number * weight, - Number * values) - { - typename Number::value_type temp[40]; - - const unsigned int points = - (fe_degree != -1 ? fe_degree : given_degree) + 1; - - AssertIndexRange(points, 40); - - 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] = 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 Number::value_type sum = 0.0; - for (unsigned int h = 0; h < points; ++h) - sum += weight[(transpose ? 1 : points) * k + - (transpose ? points : 1) * h][v] * - temp[h]; - values[i * offset + stride + j][v] = sum; - } - } - - template <int fe_degree, - unsigned int direction, - unsigned int d, - bool transpose> - static void - interpolate_3D_face(const unsigned int dof_offset, - const unsigned int given_degree, - const unsigned int v, - const Number * weight, - Number * values) - { - typename Number::value_type temp[40]; - - const unsigned int points = - (fe_degree != -1 ? fe_degree : given_degree) + 1; - - AssertIndexRange(points, 40); - - const unsigned int stride = Utilities::pow(points, 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 = 1; g < points - 1; ++g) - { - // copy result back - for (unsigned int k = 0; k < points; ++k) - temp[k] = values[dof_offset + k * stride + stride2 * g][v]; - - // perform interpolation point by point - for (unsigned int k = 0; k < points; ++k) - { - typename Number::value_type sum = 0.0; - for (unsigned int h = 0; h < points; ++h) - sum += weight[(transpose ? 1 : points) * k + - (transpose ? points : 1) * h][v] * - temp[h]; - values[dof_offset + k * stride + stride2 * g][v] = sum; - } - } - } - - template <int fe_degree, unsigned int direction, bool transpose> - static void - interpolate_3D_edge(const unsigned int p, - const unsigned int given_degree, - const unsigned int v, - const Number * weight, - Number * values) - { - typename Number::value_type temp[40]; - - const unsigned int points = - (fe_degree != -1 ? fe_degree : given_degree) + 1; - - AssertIndexRange(points, 40); - - const unsigned int stride = Utilities::pow(points, direction); - - // copy result back - for (unsigned int k = 0; k < points; ++k) - temp[k] = values[p + k * stride][v]; - - // perform interpolation point by point - for (unsigned int k = 0; k < points; ++k) - { - typename Number::value_type sum = 0.0; - for (unsigned int h = 0; h < points; ++h) - sum += weight[(transpose ? 1 : points) * k + - (transpose ? points : 1) * h][v] * - temp[h]; - values[p + k * stride][v] = sum; - } - } - - template <int fe_degree, bool transpose> - static void - run_internal(const unsigned int n_desired_components, - const FEEvaluationBaseData<dim, - typename Number::value_type, - is_face, - Number> &fe_eval, - const std::array<MatrixFreeFunctions::ConstraintKinds, - Number::size()> & 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 = - fe_eval.get_shape_info().data.front().subface_interpolation_matrices; - - const auto is_set = [](const auto a, const auto b) -> bool { - return (a & b) == b; - }; - - const auto not_set = [](const auto a, const auto b) -> bool { - return (a & b) == MatrixFreeFunctions::ConstraintKinds::unconstrained; - }; - - const unsigned int points = given_degree + 1; - - for (unsigned int c = 0; c < n_desired_components; ++c) - { - for (unsigned int v = 0; v < Number::size(); ++v) - { - const auto mask = constraint_mask[v]; - - if (mask == MatrixFreeFunctions::ConstraintKinds::unconstrained) - continue; - - if (dim == 2) // 2D: only faces - { - // direction 0: - if ((mask & MatrixFreeFunctions::ConstraintKinds::face_y) != - MatrixFreeFunctions::ConstraintKinds::unconstrained) - { - const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_x) == - MatrixFreeFunctions::ConstraintKinds::unconstrained; - - const Number *weights = - interpolation_matrices[is_subface_0].data(); - - if (is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_y)) - interpolate_2D<fe_degree, 2, transpose>( - given_degree, - v, - weights, - values); // face 2 - else - interpolate_2D<fe_degree, 3, transpose>( - given_degree, - v, - weights, - values); // face 3 - } - - // direction 1: - if ((mask & MatrixFreeFunctions::ConstraintKinds::face_x) != - MatrixFreeFunctions::ConstraintKinds::unconstrained) - { - const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_y) == - MatrixFreeFunctions::ConstraintKinds::unconstrained; - - const Number *weights = - interpolation_matrices[is_subface_0].data(); - - if (is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x)) - interpolate_2D<fe_degree, 0, transpose>( - given_degree, - v, - weights, - values); // face 0 - else - interpolate_2D<fe_degree, 1, transpose>( - given_degree, - v, - weights, - values); // face 1 - } - } - else if (dim == 3) // 3D faces and edges - { - const unsigned int p0 = 0; - const unsigned int p1 = points - 1; - const unsigned int p2 = points * points - points; - const unsigned int p3 = points * points - 1; - const unsigned int p4 = - points * points * points - points * points; - const unsigned int p5 = - points * points * points - points * points + points - 1; - const unsigned int p6 = points * points * points - points; - - const bool is_face_0 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::face_x) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x); - const bool is_face_1 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::face_x) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x); - const bool is_face_2 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::face_y) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); - const bool is_face_3 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::face_y) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); - const bool is_face_4 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::face_z) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - const bool is_face_5 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::face_z) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - - const bool is_edge_2 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_yz) && - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_y) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - const bool is_edge_3 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_yz) && - not_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_y) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - const bool is_edge_6 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_yz) && - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_y) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - const bool is_edge_7 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_yz) && - not_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_y) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - - const bool is_edge_0 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_zx) && - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - const bool is_edge_1 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_zx) && - not_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - const bool is_edge_4 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_zx) && - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - const bool is_edge_5 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_zx) && - not_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); - - const bool is_edge_8 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_xy) && - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); - const bool is_edge_9 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_xy) && - not_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); - const bool is_edge_10 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_xy) && - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); - const bool is_edge_11 = - is_set(mask, - MatrixFreeFunctions::ConstraintKinds::edge_xy) && - not_set(mask, - MatrixFreeFunctions::ConstraintKinds::type_x) && - not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); - - // direction 0: - { - const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_x) == - MatrixFreeFunctions::ConstraintKinds::unconstrained; - - const Number *weights = - interpolation_matrices[is_subface_0].data(); - - // ... faces - if (is_face_2) - interpolate_3D_face<fe_degree, 0, 1, transpose>( - p0, - given_degree, - v, - weights, - values); // face 2 - else if (is_face_3) - interpolate_3D_face<fe_degree, 0, 1, transpose>( - p2, - given_degree, - v, - weights, - values); // face 3 - if (is_face_4) - interpolate_3D_face<fe_degree, 0, 2, transpose>( - p0, - given_degree, - v, - weights, - values); // face 4 - else if (is_face_5) - interpolate_3D_face<fe_degree, 0, 2, transpose>( - p4, - given_degree, - v, - weights, - values); // face 5 - - // ... edges - if (is_face_2 || is_face_4 || is_edge_2) - interpolate_3D_edge<fe_degree, 0, transpose>( - p0, - given_degree, - v, - weights, - values); // edge 2 - if (is_face_3 || is_face_4 || is_edge_3) - interpolate_3D_edge<fe_degree, 0, transpose>( - p2, - given_degree, - v, - weights, - values); // edge 3 - if (is_face_2 || is_face_5 || is_edge_6) - interpolate_3D_edge<fe_degree, 0, transpose>( - p4, - given_degree, - v, - weights, - values); // edge 6 - if (is_face_3 || is_face_5 || is_edge_7) - interpolate_3D_edge<fe_degree, 0, transpose>( - p6, - given_degree, - v, - weights, - values); // edge 7 - } - - // direction 1: - { - const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_y) == - MatrixFreeFunctions::ConstraintKinds::unconstrained; - - const Number *weights = - interpolation_matrices[is_subface_0].data(); - - // ... faces - if (is_face_0) - interpolate_3D_face<fe_degree, 1, 0, transpose>( - p0, - given_degree, - v, - weights, - values); // face 0 - else if (is_face_1) - interpolate_3D_face<fe_degree, 1, 0, transpose>( - p1, - given_degree, - v, - weights, - values); // face 1 - if (is_face_4) - interpolate_3D_face<fe_degree, 1, 2, transpose>( - p0, - given_degree, - v, - weights, - values); // face 4 - else if (is_face_5) - interpolate_3D_face<fe_degree, 1, 2, transpose>( - p4, - given_degree, - v, - weights, - values); // face 5 - - // ... edges - if (is_face_0 || is_face_4 || is_edge_0) - interpolate_3D_edge<fe_degree, 1, transpose>( - p0, - given_degree, - v, - weights, - values); // edge 0 - if (is_face_1 || is_face_4 || is_edge_1) - interpolate_3D_edge<fe_degree, 1, transpose>( - p1, - given_degree, - v, - weights, - values); // edge 1 - if (is_face_0 || is_face_5 || is_edge_4) - interpolate_3D_edge<fe_degree, 1, transpose>( - p4, - given_degree, - v, - weights, - values); // edge 4 - if (is_face_1 || is_face_5 || is_edge_5) - interpolate_3D_edge<fe_degree, 1, transpose>( - p5, - given_degree, - v, - weights, - values); // edge 5 - } - - // direction 2: - { - const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_z) == - MatrixFreeFunctions::ConstraintKinds::unconstrained; - - const Number *weights = - interpolation_matrices[is_subface_0].data(); - - // ... faces - if (is_face_0) - interpolate_3D_face<fe_degree, 2, 0, transpose>( - p0, - given_degree, - v, - weights, - values); // face 0 - else if (is_face_1) - interpolate_3D_face<fe_degree, 2, 0, transpose>( - p1, - given_degree, - v, - weights, - values); // face 1 - if (is_face_2) - interpolate_3D_face<fe_degree, 2, 1, transpose>( - p0, - given_degree, - v, - weights, - values); // face 2 - else if (is_face_3) - interpolate_3D_face<fe_degree, 2, 1, transpose>( - p2, - given_degree, - v, - weights, - values); // face 3 - - // ... edges - if (is_face_0 || is_face_2 || is_edge_8) - interpolate_3D_edge<fe_degree, 2, transpose>( - p0, - given_degree, - v, - weights, - values); // edge 8 - if (is_face_1 || is_face_2 || is_edge_9) - interpolate_3D_edge<fe_degree, 2, transpose>( - p1, - given_degree, - v, - weights, - values); // edge 9 - if (is_face_0 || is_face_3 || is_edge_10) - interpolate_3D_edge<fe_degree, 2, transpose>( - p2, - given_degree, - v, - weights, - values); // edge 10 - if (is_face_1 || is_face_3 || is_edge_11) - interpolate_3D_edge<fe_degree, 2, transpose>( - p3, - given_degree, - v, - weights, - values); // edge 11 - } - } - else - { - Assert(false, ExcNotImplemented()); - } - } - - values += fe_eval.get_shape_info().dofs_per_component_on_cell; - } - } - }; - - } // end of namespace internal diff --git a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h new file mode 100644 index 0000000000..6b9deb6ded --- /dev/null +++ b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h @@ -0,0 +1,622 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_evaluation_kernels_hanging_nodes_h +#define dealii_matrix_free_evaluation_kernels_hanging_nodes_h + +#include <deal.II/base/config.h> + +#include <deal.II/base/vectorization.h> + +#include <deal.II/matrix_free/hanging_nodes_internal.h> + + +DEAL_II_NAMESPACE_OPEN + + +// forward declaration +template <int, typename, bool, typename> +class FEEvaluationBaseData; + + + +namespace internal +{ + template <int dim, typename Number, bool is_face> + struct FEEvaluationImplHangingNodes + { + template <int fe_degree, int n_q_points_1d> + static bool + run(const unsigned int n_desired_components, + const FEEvaluationBaseData<dim, + typename Number::value_type, + is_face, + Number> &fe_eval, + const bool transpose, + const std::array<MatrixFreeFunctions::ConstraintKinds, Number::size()> + & c_mask, + Number *values) + { + if (transpose) + run_internal<fe_degree, true>(n_desired_components, + fe_eval, + c_mask, + values); + else + run_internal<fe_degree, false>(n_desired_components, + fe_eval, + c_mask, + values); + + return false; + } + + private: + template <int fe_degree, unsigned int side, bool transpose> + static void + interpolate_2D(const unsigned int given_degree, + const unsigned int v, + const Number * weight, + Number * values) + { + typename Number::value_type temp[40]; + + const unsigned int points = + (fe_degree != -1 ? fe_degree : given_degree) + 1; + + AssertIndexRange(points, 40); + + 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] = 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 Number::value_type sum = 0.0; + for (unsigned int h = 0; h < points; ++h) + sum += weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h][v] * + temp[h]; + values[i * offset + stride + j][v] = sum; + } + } + + template <int fe_degree, + unsigned int direction, + unsigned int d, + bool transpose> + static void + interpolate_3D_face(const unsigned int dof_offset, + const unsigned int given_degree, + const unsigned int v, + const Number * weight, + Number * values) + { + typename Number::value_type temp[40]; + + const unsigned int points = + (fe_degree != -1 ? fe_degree : given_degree) + 1; + + AssertIndexRange(points, 40); + + const unsigned int stride = Utilities::pow(points, 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 = 1; g < points - 1; ++g) + { + // copy result back + for (unsigned int k = 0; k < points; ++k) + temp[k] = values[dof_offset + k * stride + stride2 * g][v]; + + // perform interpolation point by point + for (unsigned int k = 0; k < points; ++k) + { + typename Number::value_type sum = 0.0; + for (unsigned int h = 0; h < points; ++h) + sum += weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h][v] * + temp[h]; + values[dof_offset + k * stride + stride2 * g][v] = sum; + } + } + } + + template <int fe_degree, unsigned int direction, bool transpose> + static void + interpolate_3D_edge(const unsigned int p, + const unsigned int given_degree, + const unsigned int v, + const Number * weight, + Number * values) + { + typename Number::value_type temp[40]; + + const unsigned int points = + (fe_degree != -1 ? fe_degree : given_degree) + 1; + + AssertIndexRange(points, 40); + + const unsigned int stride = Utilities::pow(points, direction); + + // copy result back + for (unsigned int k = 0; k < points; ++k) + temp[k] = values[p + k * stride][v]; + + // perform interpolation point by point + for (unsigned int k = 0; k < points; ++k) + { + typename Number::value_type sum = 0.0; + for (unsigned int h = 0; h < points; ++h) + sum += weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h][v] * + temp[h]; + values[p + k * stride][v] = sum; + } + } + + template <int fe_degree, bool transpose> + static void + run_internal(const unsigned int n_desired_components, + const FEEvaluationBaseData<dim, + typename Number::value_type, + is_face, + Number> &fe_eval, + const std::array<MatrixFreeFunctions::ConstraintKinds, + Number::size()> & 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 = + fe_eval.get_shape_info().data.front().subface_interpolation_matrices; + + const auto is_set = [](const auto a, const auto b) -> bool { + return (a & b) == b; + }; + + const auto not_set = [](const auto a, const auto b) -> bool { + return (a & b) == MatrixFreeFunctions::ConstraintKinds::unconstrained; + }; + + const unsigned int points = given_degree + 1; + + for (unsigned int c = 0; c < n_desired_components; ++c) + { + for (unsigned int v = 0; v < Number::size(); ++v) + { + const auto mask = constraint_mask[v]; + + if (mask == MatrixFreeFunctions::ConstraintKinds::unconstrained) + continue; + + if (dim == 2) // 2D: only faces + { + // direction 0: + if ((mask & MatrixFreeFunctions::ConstraintKinds::face_y) != + MatrixFreeFunctions::ConstraintKinds::unconstrained) + { + const bool is_subface_0 = + (mask & MatrixFreeFunctions::ConstraintKinds::type_x) == + MatrixFreeFunctions::ConstraintKinds::unconstrained; + + const Number *weights = + interpolation_matrices[is_subface_0].data(); + + if (is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_y)) + interpolate_2D<fe_degree, 2, transpose>( + given_degree, + v, + weights, + values); // face 2 + else + interpolate_2D<fe_degree, 3, transpose>( + given_degree, + v, + weights, + values); // face 3 + } + + // direction 1: + if ((mask & MatrixFreeFunctions::ConstraintKinds::face_x) != + MatrixFreeFunctions::ConstraintKinds::unconstrained) + { + const bool is_subface_0 = + (mask & MatrixFreeFunctions::ConstraintKinds::type_y) == + MatrixFreeFunctions::ConstraintKinds::unconstrained; + + const Number *weights = + interpolation_matrices[is_subface_0].data(); + + if (is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x)) + interpolate_2D<fe_degree, 0, transpose>( + given_degree, + v, + weights, + values); // face 0 + else + interpolate_2D<fe_degree, 1, transpose>( + given_degree, + v, + weights, + values); // face 1 + } + } + else if (dim == 3) // 3D faces and edges + { + const unsigned int p0 = 0; + const unsigned int p1 = points - 1; + const unsigned int p2 = points * points - points; + const unsigned int p3 = points * points - 1; + const unsigned int p4 = + points * points * points - points * points; + const unsigned int p5 = + points * points * points - points * points + points - 1; + const unsigned int p6 = points * points * points - points; + + const bool is_face_0 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::face_x) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x); + const bool is_face_1 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::face_x) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x); + const bool is_face_2 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::face_y) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); + const bool is_face_3 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::face_y) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); + const bool is_face_4 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::face_z) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + const bool is_face_5 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::face_z) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + + const bool is_edge_2 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_yz) && + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_y) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + const bool is_edge_3 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_yz) && + not_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_y) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + const bool is_edge_6 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_yz) && + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_y) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + const bool is_edge_7 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_yz) && + not_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_y) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + + const bool is_edge_0 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_zx) && + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + const bool is_edge_1 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_zx) && + not_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + const bool is_edge_4 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_zx) && + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + const bool is_edge_5 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_zx) && + not_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z); + + const bool is_edge_8 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_xy) && + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); + const bool is_edge_9 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_xy) && + not_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); + const bool is_edge_10 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_xy) && + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); + const bool is_edge_11 = + is_set(mask, + MatrixFreeFunctions::ConstraintKinds::edge_xy) && + not_set(mask, + MatrixFreeFunctions::ConstraintKinds::type_x) && + not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y); + + // direction 0: + { + const bool is_subface_0 = + (mask & MatrixFreeFunctions::ConstraintKinds::type_x) == + MatrixFreeFunctions::ConstraintKinds::unconstrained; + + const Number *weights = + interpolation_matrices[is_subface_0].data(); + + // ... faces + if (is_face_2) + interpolate_3D_face<fe_degree, 0, 1, transpose>( + p0, + given_degree, + v, + weights, + values); // face 2 + else if (is_face_3) + interpolate_3D_face<fe_degree, 0, 1, transpose>( + p2, + given_degree, + v, + weights, + values); // face 3 + if (is_face_4) + interpolate_3D_face<fe_degree, 0, 2, transpose>( + p0, + given_degree, + v, + weights, + values); // face 4 + else if (is_face_5) + interpolate_3D_face<fe_degree, 0, 2, transpose>( + p4, + given_degree, + v, + weights, + values); // face 5 + + // ... edges + if (is_face_2 || is_face_4 || is_edge_2) + interpolate_3D_edge<fe_degree, 0, transpose>( + p0, + given_degree, + v, + weights, + values); // edge 2 + if (is_face_3 || is_face_4 || is_edge_3) + interpolate_3D_edge<fe_degree, 0, transpose>( + p2, + given_degree, + v, + weights, + values); // edge 3 + if (is_face_2 || is_face_5 || is_edge_6) + interpolate_3D_edge<fe_degree, 0, transpose>( + p4, + given_degree, + v, + weights, + values); // edge 6 + if (is_face_3 || is_face_5 || is_edge_7) + interpolate_3D_edge<fe_degree, 0, transpose>( + p6, + given_degree, + v, + weights, + values); // edge 7 + } + + // direction 1: + { + const bool is_subface_0 = + (mask & MatrixFreeFunctions::ConstraintKinds::type_y) == + MatrixFreeFunctions::ConstraintKinds::unconstrained; + + const Number *weights = + interpolation_matrices[is_subface_0].data(); + + // ... faces + if (is_face_0) + interpolate_3D_face<fe_degree, 1, 0, transpose>( + p0, + given_degree, + v, + weights, + values); // face 0 + else if (is_face_1) + interpolate_3D_face<fe_degree, 1, 0, transpose>( + p1, + given_degree, + v, + weights, + values); // face 1 + if (is_face_4) + interpolate_3D_face<fe_degree, 1, 2, transpose>( + p0, + given_degree, + v, + weights, + values); // face 4 + else if (is_face_5) + interpolate_3D_face<fe_degree, 1, 2, transpose>( + p4, + given_degree, + v, + weights, + values); // face 5 + + // ... edges + if (is_face_0 || is_face_4 || is_edge_0) + interpolate_3D_edge<fe_degree, 1, transpose>( + p0, + given_degree, + v, + weights, + values); // edge 0 + if (is_face_1 || is_face_4 || is_edge_1) + interpolate_3D_edge<fe_degree, 1, transpose>( + p1, + given_degree, + v, + weights, + values); // edge 1 + if (is_face_0 || is_face_5 || is_edge_4) + interpolate_3D_edge<fe_degree, 1, transpose>( + p4, + given_degree, + v, + weights, + values); // edge 4 + if (is_face_1 || is_face_5 || is_edge_5) + interpolate_3D_edge<fe_degree, 1, transpose>( + p5, + given_degree, + v, + weights, + values); // edge 5 + } + + // direction 2: + { + const bool is_subface_0 = + (mask & MatrixFreeFunctions::ConstraintKinds::type_z) == + MatrixFreeFunctions::ConstraintKinds::unconstrained; + + const Number *weights = + interpolation_matrices[is_subface_0].data(); + + // ... faces + if (is_face_0) + interpolate_3D_face<fe_degree, 2, 0, transpose>( + p0, + given_degree, + v, + weights, + values); // face 0 + else if (is_face_1) + interpolate_3D_face<fe_degree, 2, 0, transpose>( + p1, + given_degree, + v, + weights, + values); // face 1 + if (is_face_2) + interpolate_3D_face<fe_degree, 2, 1, transpose>( + p0, + given_degree, + v, + weights, + values); // face 2 + else if (is_face_3) + interpolate_3D_face<fe_degree, 2, 1, transpose>( + p2, + given_degree, + v, + weights, + values); // face 3 + + // ... edges + if (is_face_0 || is_face_2 || is_edge_8) + interpolate_3D_edge<fe_degree, 2, transpose>( + p0, + given_degree, + v, + weights, + values); // edge 8 + if (is_face_1 || is_face_2 || is_edge_9) + interpolate_3D_edge<fe_degree, 2, transpose>( + p1, + given_degree, + v, + weights, + values); // edge 9 + if (is_face_0 || is_face_3 || is_edge_10) + interpolate_3D_edge<fe_degree, 2, transpose>( + p2, + given_degree, + v, + weights, + values); // edge 10 + if (is_face_1 || is_face_3 || is_edge_11) + interpolate_3D_edge<fe_degree, 2, transpose>( + p3, + given_degree, + v, + weights, + values); // edge 11 + } + } + else + { + Assert(false, ExcNotImplemented()); + } + } + + values += fe_eval.get_shape_info().dofs_per_component_on_cell; + } + } + }; + + +} // end of namespace internal + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h index 62fc91bcd0..1c9d8b61ef 100644 --- a/include/deal.II/matrix_free/evaluation_selector.h +++ b/include/deal.II/matrix_free/evaluation_selector.h @@ -19,7 +19,9 @@ #include <deal.II/base/config.h> +#include <deal.II/matrix_free/evaluation_flags.h> #include <deal.II/matrix_free/evaluation_kernels.h> +#include <deal.II/matrix_free/shape_info.h> DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index f242cb5d23..b246ee4e44 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -23,46 +23,13 @@ #include <deal.II/matrix_free/evaluation_kernels.h> #include <deal.II/matrix_free/evaluation_selector.h> #include <deal.II/matrix_free/evaluation_template_factory.h> +#include <deal.II/matrix_free/evaluation_template_factory_internal.h> #include <deal.II/matrix_free/fe_evaluation.h> -#ifndef FE_EVAL_FACTORY_DEGREE_MAX -# define FE_EVAL_FACTORY_DEGREE_MAX 6 -#endif - DEAL_II_NAMESPACE_OPEN namespace internal { - template <int degree, typename EvaluatorType, typename... Args> - bool - instantiation_helper_run(const unsigned int given_degree, - const unsigned int n_q_points_1d, - Args &...args) - { - if (given_degree == degree) - { - if (n_q_points_1d == degree + 1) - return EvaluatorType::template run<degree, degree + 1>(args...); - else if (n_q_points_1d == degree + 2) - return EvaluatorType::template run<degree, degree + 2>(args...); - else if (n_q_points_1d == degree) - return EvaluatorType::template run<degree, degree>(args...); - else if (n_q_points_1d == (3 * degree) / 2 + 1) - return EvaluatorType::template run<degree, (3 * degree) / 2 + 1>( - args...); - else - // slow path - return EvaluatorType::template run<-1, 0>(args...); - } - else if (degree < FE_EVAL_FACTORY_DEGREE_MAX) - return instantiation_helper_run< - (degree < FE_EVAL_FACTORY_DEGREE_MAX ? degree + 1 : degree), - EvaluatorType>(given_degree, n_q_points_1d, args...); - else - // slow path - return EvaluatorType::template run<-1, 0>(args...); - } - struct FastEvaluationSupported { template <int fe_degree, int n_q_points_1d> @@ -420,57 +387,6 @@ namespace internal fe_degree, n_q_points_1d, n_components, fe_eval, in_array, out_array); } - - - template <int dim, typename Number, typename VectorizedArrayType> - void - FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply( - const unsigned int n_components, - const unsigned int fe_degree, - const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> - & fe_eval, - const bool transpose, - const std::array<MatrixFreeFunctions::ConstraintKinds, - VectorizedArrayType::size()> &c_mask, - VectorizedArrayType * values) - { - instantiation_helper_run< - 1, - FEEvaluationImplHangingNodes<dim, VectorizedArrayType, false>>( - fe_degree, - fe_degree + 1, - n_components, - fe_eval, - transpose, - c_mask, - values); - } - - - - template <int dim, typename Number, typename VectorizedArrayType> - void - FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply( - const unsigned int n_components, - const unsigned int fe_degree, - const FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &fe_eval, - const bool transpose, - const std::array<MatrixFreeFunctions::ConstraintKinds, - VectorizedArrayType::size()> &c_mask, - VectorizedArrayType * values) - { - instantiation_helper_run< - 1, - FEEvaluationImplHangingNodes<dim, VectorizedArrayType, true>>( - fe_degree, - fe_degree + 1, - n_components, - fe_eval, - transpose, - c_mask, - values); - } - } // end of namespace internal DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h new file mode 100644 index 0000000000..cf03c7c2d5 --- /dev/null +++ b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_evaluation_template_factory_hanging_nodes_templates_h +#define dealii_matrix_free_evaluation_template_factory_hanging_nodes_templates_h + + +#include <deal.II/base/config.h> + +#include <deal.II/matrix_free/evaluation_kernels_hanging_nodes.h> +#include <deal.II/matrix_free/evaluation_template_factory.h> +#include <deal.II/matrix_free/evaluation_template_factory_internal.h> +#include <deal.II/matrix_free/fe_evaluation.h> + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + template <int dim, typename Number, typename VectorizedArrayType> + void + FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply( + const unsigned int n_components, + const unsigned int fe_degree, + const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> + & fe_eval, + const bool transpose, + const std::array<MatrixFreeFunctions::ConstraintKinds, + VectorizedArrayType::size()> &c_mask, + VectorizedArrayType * values) + { + instantiation_helper_degree_run< + 1, + FEEvaluationImplHangingNodes<dim, VectorizedArrayType, false>>( + fe_degree, n_components, fe_eval, transpose, c_mask, values); + } + + + + template <int dim, typename Number, typename VectorizedArrayType> + void + FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply( + const unsigned int n_components, + const unsigned int fe_degree, + const FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &fe_eval, + const bool transpose, + const std::array<MatrixFreeFunctions::ConstraintKinds, + VectorizedArrayType::size()> &c_mask, + VectorizedArrayType * values) + { + instantiation_helper_degree_run< + 1, + FEEvaluationImplHangingNodes<dim, VectorizedArrayType, true>>( + fe_degree, n_components, fe_eval, transpose, c_mask, values); + } +} // end of namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/evaluation_template_factory_internal.h b/include/deal.II/matrix_free/evaluation_template_factory_internal.h new file mode 100644 index 0000000000..94f66caead --- /dev/null +++ b/include/deal.II/matrix_free/evaluation_template_factory_internal.h @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_evaluation_template_factory_internal_h +#define dealii_matrix_free_evaluation_template_factory_internal_h + + +#include <deal.II/base/config.h> + +#ifndef FE_EVAL_FACTORY_DEGREE_MAX +# define FE_EVAL_FACTORY_DEGREE_MAX 6 +#endif + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + template <int degree, typename EvaluatorType, typename... Args> + bool + instantiation_helper_run(const unsigned int given_degree, + const unsigned int n_q_points_1d, + Args &...args) + { + if (given_degree == degree) + { + if (n_q_points_1d == degree + 1) + return EvaluatorType::template run<degree, degree + 1>(args...); + else if (n_q_points_1d == degree + 2) + return EvaluatorType::template run<degree, degree + 2>(args...); + else if (n_q_points_1d == degree) + return EvaluatorType::template run<degree, degree>(args...); + else if (n_q_points_1d == (3 * degree) / 2 + 1) + return EvaluatorType::template run<degree, (3 * degree) / 2 + 1>( + args...); + else + // slow path + return EvaluatorType::template run<-1, 0>(args...); + } + else if (degree < FE_EVAL_FACTORY_DEGREE_MAX) + return instantiation_helper_run< + (degree < FE_EVAL_FACTORY_DEGREE_MAX ? degree + 1 : degree), + EvaluatorType>(given_degree, n_q_points_1d, args...); + else + // slow path + return EvaluatorType::template run<-1, 0>(args...); + } + + template <int degree, typename EvaluatorType, typename... Args> + bool + instantiation_helper_degree_run(const unsigned int given_degree, + Args &...args) + { + if (given_degree == degree) + { + return EvaluatorType::template run<degree, degree + 1>(args...); + } + else if (degree < FE_EVAL_FACTORY_DEGREE_MAX) + return instantiation_helper_degree_run< + (degree < FE_EVAL_FACTORY_DEGREE_MAX ? degree + 1 : degree), + EvaluatorType>(given_degree, args...); + else + // slow path + return EvaluatorType::template run<-1, 0>(args...); + } + +} // end of namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/matrix_free/CMakeLists.txt b/source/matrix_free/CMakeLists.txt index bb85bc75d7..716b9b8511 100644 --- a/source/matrix_free/CMakeLists.txt +++ b/source/matrix_free/CMakeLists.txt @@ -23,6 +23,12 @@ SET(_src evaluation_template_factory_inst4.cc evaluation_template_factory_inst5.cc evaluation_template_factory_inst6.cc + evaluation_template_factory_hanging_nodes.cc + evaluation_template_factory_hanging_nodes_inst2.cc + evaluation_template_factory_hanging_nodes_inst3.cc + evaluation_template_factory_hanging_nodes_inst4.cc + evaluation_template_factory_hanging_nodes_inst5.cc + evaluation_template_factory_hanging_nodes_inst6.cc fe_point_evaluation.cc mapping_info.cc mapping_info_inst2.cc @@ -35,6 +41,7 @@ SET(_src SET(_inst evaluation_template_factory.inst.in + evaluation_template_factory_hanging_nodes.inst.in fe_point_evaluation.inst.in mapping_info.inst.in matrix_free.inst.in diff --git a/source/matrix_free/evaluation_template_factory.inst.in b/source/matrix_free/evaluation_template_factory.inst.in index f3e05fbf09..468f6f3724 100644 --- a/source/matrix_free/evaluation_template_factory.inst.in +++ b/source/matrix_free/evaluation_template_factory.inst.in @@ -32,10 +32,4 @@ for (deal_II_dimension : DIMENSIONS; deal_II_dimension, deal_II_scalar_vectorized::value_type, deal_II_scalar_vectorized>; - - // inverse mass - template struct dealii::internal::FEEvaluationHangingNodesFactory< - deal_II_dimension, - deal_II_scalar_vectorized::value_type, - deal_II_scalar_vectorized>; } diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes.cc new file mode 100644 index 0000000000..0d7a73694f --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_hanging_nodes.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include <deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h> + +DEAL_II_NAMESPACE_OPEN + +#define SPLIT_INSTANTIATIONS_COUNT 6 +#ifndef SPLIT_INSTANTIATIONS_INDEX +# define SPLIT_INSTANTIATIONS_INDEX 0 +#endif + +#include "evaluation_template_factory_hanging_nodes.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes.inst.in b/source/matrix_free/evaluation_template_factory_hanging_nodes.inst.in new file mode 100644 index 0000000000..0ed6662cae --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_hanging_nodes.inst.in @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS; + deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) + { + // hanging nodes + template struct dealii::internal::FEEvaluationHangingNodesFactory< + deal_II_dimension, + deal_II_scalar_vectorized::value_type, + deal_II_scalar_vectorized>; + } diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst2.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst2.cc new file mode 100644 index 0000000000..0a69454668 --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst2.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 1 +#include "evaluation_template_factory_hanging_nodes.cc" diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst3.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst3.cc new file mode 100644 index 0000000000..4e1742578b --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst3.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 2 +#include "evaluation_template_factory_hanging_nodes.cc" diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst4.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst4.cc new file mode 100644 index 0000000000..cab92f149f --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst4.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 3 +#include "evaluation_template_factory_hanging_nodes.cc" diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst5.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst5.cc new file mode 100644 index 0000000000..ea2904f3f2 --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst5.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 4 +#include "evaluation_template_factory_hanging_nodes.cc" diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst6.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst6.cc new file mode 100644 index 0000000000..e6d89c69ce --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst6.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 5 +#include "evaluation_template_factory_hanging_nodes.cc" diff --git a/tests/matrix_free/hanging_node_kernels_01.cc b/tests/matrix_free/hanging_node_kernels_01.cc index 4ee4fe5097..d93b708336 100644 --- a/tests/matrix_free/hanging_node_kernels_01.cc +++ b/tests/matrix_free/hanging_node_kernels_01.cc @@ -24,6 +24,7 @@ #include <deal.II/grid/grid_generator.h> +#include <deal.II/matrix_free/evaluation_kernels_hanging_nodes.h> #include <deal.II/matrix_free/fe_evaluation.h> #include <deal.II/matrix_free/matrix_free.h>