#include <deal.II/base/vectorization.h>
#include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/evaluation_kernels_common.h>
#include <deal.II/matrix_free/fe_evaluation_data.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>
(type == MatrixFreeFunctions::truncated_tensor) ?
Utilities::pow(shape_data.front().fe_degree + 1, dim) :
fe_eval.get_shape_info().dofs_per_component_on_cell;
- const Number *values_dofs = values_dofs_actual;
+ const Number *values_dofs =
+ (type == MatrixFreeFunctions::truncated_tensor) ?
+ temp1 + 2 * (std::max<std::size_t>(
+ fe_eval.get_shape_info().dofs_per_component_on_cell,
+ n_q_points)) :
+ values_dofs_actual;
+
if (type == MatrixFreeFunctions::truncated_tensor)
- {
- const std::size_t n_dofs_per_comp =
- fe_eval.get_shape_info().dofs_per_component_on_cell;
- Number *values_dofs_tmp =
- temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points));
- const int degree =
- fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
- for (unsigned int c = 0; c < n_components; ++c)
- for (int i = 0, count_p = 0, count_q = 0;
- i < (dim > 2 ? degree + 1 : 1);
- ++i)
- {
- for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
- {
- for (int k = 0; k < degree + 1 - j - i;
- ++k, ++count_p, ++count_q)
- values_dofs_tmp[c * dofs_per_comp + count_q] =
- values_dofs_actual[c * n_dofs_per_comp + count_p];
- for (int k = degree + 1 - j - i; k < degree + 1;
- ++k, ++count_q)
- values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
- }
- for (int j = degree + 1 - i; j < degree + 1; ++j)
- for (int k = 0; k < degree + 1; ++k, ++count_q)
- values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
- }
- values_dofs = values_dofs_tmp;
- }
+ embed_truncated_into_full_tensor_product<dim, fe_degree>(
+ n_components,
+ const_cast<Number *>(values_dofs),
+ values_dofs_actual,
+ fe_eval);
Number *values_quad = fe_eval.begin_values();
Number *gradients_quad = fe_eval.begin_gradients();
(type == MatrixFreeFunctions::truncated_tensor) ?
Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
fe_eval.get_shape_info().dofs_per_component_on_cell;
+
// expand dof_values to tensor product for truncated tensor products
Number *values_dofs =
(type == MatrixFreeFunctions::truncated_tensor) ?
}
if (type == MatrixFreeFunctions::truncated_tensor)
- {
- const std::size_t n_dofs_per_comp =
- fe_eval.get_shape_info().dofs_per_component_on_cell;
- values_dofs -= dofs_per_comp * n_components;
- const int degree =
- fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
- for (unsigned int c = 0; c < n_components; ++c)
- for (int i = 0, count_p = 0, count_q = 0;
- i < (dim > 2 ? degree + 1 : 1);
- ++i)
- {
- for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
- {
- for (int k = 0; k < degree + 1 - j - i;
- ++k, ++count_p, ++count_q)
- values_dofs_actual[c * n_dofs_per_comp + count_p] =
- values_dofs[c * dofs_per_comp + count_q];
- count_q += j + i;
- }
- count_q += i * (degree + 1);
- }
- }
+ truncate_tensor_product_to_complete_degrees<dim, fe_degree>(
+ n_components,
+ values_dofs_actual,
+ values_dofs - dofs_per_comp * n_components,
+ fe_eval);
}
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_evaluation_kernels_common_h
+#define dealii_matrix_free_evaluation_kernels_common_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/matrix_free/fe_evaluation_data.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+ template <int dim, int fe_degree, typename Number, bool is_face>
+ inline void
+ embed_truncated_into_full_tensor_product(
+ const unsigned int n_components,
+ Number *values_dofs,
+ const Number *values_dofs_actual,
+ FEEvaluationData<dim, Number, is_face> &fe_eval)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const std::size_t dofs_per_comp =
+ Utilities::pow(shape_data.fe_degree + 1, dim);
+ const std::size_t n_dofs_per_comp = shape_info.dofs_per_component_on_cell;
+ const int degree = fe_degree != -1 ? fe_degree : shape_data.fe_degree;
+
+ for (unsigned int c = 0; c < n_components; ++c)
+ for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1);
+ ++i)
+ {
+ for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
+ {
+ for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q)
+ values_dofs[c * dofs_per_comp + count_q] =
+ values_dofs_actual[c * n_dofs_per_comp + count_p];
+ for (int k = degree + 1 - j - i; k < degree + 1; ++k, ++count_q)
+ values_dofs[c * dofs_per_comp + count_q] = Number();
+ }
+ for (int j = degree + 1 - i; j < degree + 1; ++j)
+ for (int k = 0; k < degree + 1; ++k, ++count_q)
+ values_dofs[c * dofs_per_comp + count_q] = Number();
+ }
+ }
+
+
+
+ template <int dim, int fe_degree, typename Number, bool is_face>
+ inline void
+ truncate_tensor_product_to_complete_degrees(
+ const unsigned int n_components,
+ Number *values_dofs_actual,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, is_face> &fe_eval)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp =
+ Utilities::fixed_power<dim>(shape_data.fe_degree + 1);
+ const std::size_t n_dofs_per_comp = shape_info.dofs_per_component_on_cell;
+ const int degree = fe_degree != -1 ? fe_degree : shape_data.fe_degree;
+ for (unsigned int c = 0; c < n_components; ++c)
+ for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1);
+ ++i)
+ {
+ for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
+ {
+ for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q)
+ values_dofs_actual[c * n_dofs_per_comp + count_p] =
+ values_dofs[c * dofs_per_comp + count_q];
+ count_q += j + i;
+ }
+ count_q += i * (degree + 1);
+ }
+ }
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <deal.II/matrix_free/dof_info.h>
#include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/evaluation_kernels_common.h>
#include <deal.II/matrix_free/fe_evaluation_data.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>
interpolate_raviart_thomas<do_evaluate, add_into_output>(
n_components, input, output, flags, face_no, shape_info);
else
- interpolate_generic<do_evaluate, add_into_output>(
- n_components,
- input,
- output,
- flags,
- face_no,
- shape_info.data.front().fe_degree + 1,
- shape_info.data.front().shape_data_on_face,
- shape_info.dofs_per_component_on_cell,
- 3 * shape_info.dofs_per_component_on_face);
+ {
+ const unsigned int fe_degree_ = shape_info.data.front().fe_degree;
+
+ interpolate_generic<do_evaluate, add_into_output>(
+ n_components,
+ input,
+ output,
+ flags,
+ face_no,
+ fe_degree_ + 1,
+ shape_info.data.front().shape_data_on_face,
+ Utilities::pow(fe_degree_ + 1, dim),
+ 3 * Utilities::pow(fe_degree_ + 1, dim - 1));
+ }
}
/**
static bool
evaluate_tensor(const unsigned int n_components,
const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number *values_dofs,
+ const Number *values_dofs_actual,
FEEvaluationData<dim, Number, true> &fe_eval)
{
const auto &shape_info = fe_eval.get_shape_info();
// Note: we always keep storage of values, 1st and 2nd derivatives in an
// array, so reserve space for all three here
- Number *temp = fe_eval.get_scratch_data().begin();
- Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+ Number *temp1 = fe_eval.get_scratch_data().begin();
+ Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
+
+ const Number *values_dofs =
+ (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
+ temp2 +
+ 2 * (std::max(fe_eval.get_shape_info().dofs_per_component_on_cell,
+ shape_info.n_q_points)) :
+ values_dofs_actual;
+
+ if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
+ embed_truncated_into_full_tensor_product<dim, fe_degree>(
+ n_components,
+ const_cast<Number *>(values_dofs),
+ values_dofs_actual,
+ fe_eval);
bool use_vectorization = true;
if (fe_eval.get_dof_access_index() ==
values_dofs,
fe_eval,
use_vectorization,
- temp,
- scratch_data);
+ temp1,
+ temp2);
evaluate_in_face<fe_degree, n_q_points_1d>(
- n_components, evaluation_flag, fe_eval, temp, scratch_data);
+ n_components, evaluation_flag, fe_eval, temp1, temp2);
if (dim == 3)
adjust_quadrature_for_face_orientation(
- n_components, evaluation_flag, fe_eval, use_vectorization, temp);
+ n_components, evaluation_flag, fe_eval, use_vectorization, temp1);
return false;
}
static bool
integrate_tensor(const unsigned int n_components,
const EvaluationFlags::EvaluationFlags integration_flag,
- Number *values_dofs,
+ Number *values_dofs_actual,
FEEvaluationData<dim, Number, true> &fe_eval,
const bool sum_into_values)
{
Utilities::pow(fe_degree + 1, dim - 1) :
Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
- Number *temp = fe_eval.get_scratch_data().begin();
- Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+ Number *temp1 = fe_eval.get_scratch_data().begin();
+ Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
+
+ // expand dof_values to tensor product for truncated tensor products
+ Number *values_dofs =
+ (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
+ temp2 + 2 * (std::max<std::size_t>(
+ fe_eval.get_shape_info().dofs_per_component_on_cell,
+ fe_eval.get_shape_info().n_q_points)) :
+ values_dofs_actual;
bool use_vectorization = true;
if (dim == 3)
adjust_quadrature_for_face_orientation(
- n_components, integration_flag, fe_eval, use_vectorization, temp);
+ n_components, integration_flag, fe_eval, use_vectorization, temp1);
integrate_in_face<fe_degree, n_q_points_1d>(
- n_components, integration_flag, fe_eval, temp, scratch_data);
+ n_components, integration_flag, fe_eval, temp1, temp2);
collect_from_face<fe_degree>(n_components,
integration_flag,
values_dofs,
fe_eval,
use_vectorization,
- temp,
- scratch_data,
+ temp1,
+ temp2,
sum_into_values);
+
+ if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
+ truncate_tensor_product_to_complete_degrees<dim, fe_degree>(
+ n_components, values_dofs_actual, values_dofs, fe_eval);
+
return false;
}
const unsigned int size_scratch_data =
std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components *
- 3 +
+ 4 +
2 * n_quadrature_points;
const unsigned int size_data_arrays =
n_components * dofs_per_component +
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// Similar to matrix_vector_faces_36 but with FE_DGP.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/fe/fe_dgp.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_faces_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(5 - dim);
+
+ FE_DGP<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ do_test<dim, fe_degree, fe_degree + 1, double>(dof, constraints, true);
+}
--- /dev/null
+
+DEAL:2d::Testing FE_DGP<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::Norm of difference parallel: 0
+DEAL:2d::
+DEAL:2d::Testing FE_DGP<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::Norm of difference parallel: 0
+DEAL:2d::
+DEAL:3d::Testing FE_DGP<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::Norm of difference parallel: 0
+DEAL:3d::
+DEAL:3d::Testing FE_DGP<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::Norm of difference parallel: 0
+DEAL:3d::