From 29c7b3f3178e23c3d6a2ed8060ee3002b622560c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 15 Nov 2020 17:57:46 +0100 Subject: [PATCH] Make FEEvaluationImpl work for ElementType::tensor_none --- .../deal.II/matrix_free/evaluation_kernels.h | 180 ++++++++++++++++++ include/deal.II/matrix_free/shape_info.h | 7 +- 2 files changed, 186 insertions(+), 1 deletion(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index eb9c18af22..036b628f59 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -137,6 +137,40 @@ namespace internal + /** + * Specialization for MatrixFreeFunctions::tensor_none, which cannot use the + * sum-factorization kernels. + */ + template + struct FEEvaluationImpl + { + static void + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const Number * values_dofs_actual, + Number * values_quad, + Number * gradients_quad, + Number * hessians_quad, + Number * scratch_data); + + static void + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const MatrixFreeFunctions::ShapeInfo &shape_info, + Number * values_dofs_actual, + Number * values_quad, + Number * gradients_quad, + Number * scratch_data, + const bool add_into_values_array); + }; + + + template + inline void + FEEvaluationImpl< + MatrixFreeFunctions::tensor_none, + dim, + fe_degree, + n_q_points_1d, + Number>::evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const Number *values_dofs_actual, + Number * values_quad, + Number * gradients_quad, + Number * hessians_quad, + Number * scratch_data) + { + (void)scratch_data; + + const unsigned int n_dofs = shape_info.dofs_per_component_on_cell; + const unsigned int n_q_points = shape_info.n_q_points; + + if (evaluation_flag & EvaluationFlags::values) + { + const auto shape_values = shape_info.data.front().shape_values.data(); + for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int q = 0; q < n_q_points; ++q) + { + values_quad[n_q_points * c + q] = 0.0; + for (unsigned int i = 0; i < n_dofs; ++i) + values_quad[n_q_points * c + q] += + shape_values[i * n_q_points + q] * + values_dofs_actual[n_dofs * c + i]; + } + } + + if (evaluation_flag & EvaluationFlags::gradients) + { + const auto shape_gradients = + shape_info.data.front().shape_gradients.data(); + for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int q = 0; q < n_q_points; ++q) + { + gradients_quad[n_q_points * dim * c + n_q_points * d + q] = 0.0; + for (unsigned int i = 0; i < n_dofs; ++i) + gradients_quad[n_q_points * dim * c + n_q_points * d + q] += + shape_gradients[d * n_dofs * n_q_points + i * n_q_points + + q] * + values_dofs_actual[n_dofs * dim * c + i]; + } + } + + if (evaluation_flag & EvaluationFlags::hessians) + { + Assert(false, ExcNotImplemented()); + (void)hessians_quad; + } + } + + + + template + inline void + FEEvaluationImpl< + MatrixFreeFunctions::tensor_none, + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const MatrixFreeFunctions::ShapeInfo &shape_info, + Number * values_dofs_actual, + Number * values_quad, + Number * gradients_quad, + Number * scratch_data, + const bool add_into_values_array) + { + (void)scratch_data; + + const unsigned int n_dofs = shape_info.dofs_per_component_on_cell; + const unsigned int n_q_points = shape_info.n_q_points; + + if (add_into_values_array == false) + for (unsigned int i = 0; i < n_components * n_dofs; ++i) + values_dofs_actual[i] = 0.0; + + if (integration_flag & EvaluationFlags::values) + { + const auto shape_values = shape_info.data.front().shape_values.data(); + for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int i = 0; i < n_dofs; ++i) + values_dofs_actual[n_dofs * c + i] += + shape_values[i * n_q_points + q] * + values_quad[n_q_points * c + q]; + } + + if (integration_flag & EvaluationFlags::gradients) + { + const auto shape_gradients = + shape_info.data.front().shape_gradients.data(); + for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int i = 0; i < n_dofs; ++i) + values_dofs_actual[n_dofs * dim * c + i] += + shape_gradients[d * n_dofs * n_q_points + i * n_q_points + + q] * + gradients_quad[n_q_points * dim * c + n_q_points * d + q]; + } + } + + + /** * This struct implements the change between two different bases. This is an * ingredient in the FEEvaluationImplTransformToCollocation class where we @@ -1499,6 +1647,22 @@ namespace internal hessians_quad, scratch_data); } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::tensor_none) + { + internal::FEEvaluationImpl::evaluate(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); + } else { internal::FEEvaluationImpl< @@ -1645,6 +1809,22 @@ namespace internal scratch_data, sum_into_values_array); } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::tensor_none) + { + internal::FEEvaluationImpl::integrate(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); + } else { internal::FEEvaluationImpl< diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index bb896a9c54..12f3b6eea1 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -91,7 +91,12 @@ namespace internal * of the unit interval 0.5 that additionally add a constant shape * function according to FE_Q_DG0. */ - tensor_symmetric_plus_dg0 = 5 + tensor_symmetric_plus_dg0 = 5, + + /** + * Shape functions without an tensor product properties. + */ + tensor_none = 6 }; -- 2.39.5