From 113a3148736057331e93425b822827379aee5004 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 1 Sep 2020 12:48:10 +0200 Subject: [PATCH] Work on interpolation from cell to face quadrature points --- .../deal.II/matrix_free/evaluation_kernels.h | 69 +++++++- include/deal.II/matrix_free/shape_info.h | 12 ++ .../matrix_free/shape_info.templates.h | 25 ++- .../matrix_free/interpolate_quadrature_01.cc | 157 ++++++++++++++++++ .../interpolate_quadrature_01.output | 2 + 5 files changed, 255 insertions(+), 10 deletions(-) create mode 100644 tests/matrix_free/interpolate_quadrature_01.cc create mode 100644 tests/matrix_free/interpolate_quadrature_01.output diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 4451027de4..aaad644911 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -1686,24 +1686,77 @@ namespace internal Number * output, const bool do_gradients, const unsigned int face_no) + { + Assert(static_cast(fe_degree) == + data.data.front().fe_degree || + fe_degree == -1, + ExcInternalError()); + + interpolate_generic( + input, + output, + do_gradients, + face_no, + data.data.front().fe_degree + 1, + data.data.front().shape_data_on_face, + data.dofs_per_component_on_cell, + 2 * data.dofs_per_component_on_face); + } + + /** + * Interpolate the values on the cell quadrature points onto a face. + */ + template + static void + interpolate_quadrature(const MatrixFreeFunctions::ShapeInfo &data, + const Number * input, + Number * output, + const bool do_gradients, + const unsigned int face_no) + { + Assert(static_cast(fe_degree + 1) == + data.data.front().quadrature.size() || + fe_degree == -1, + ExcInternalError()); + + interpolate_generic( + input, + output, + do_gradients, + face_no, + data.data.front().quadrature.size(), + data.data.front().quadrature_data_on_face, + data.n_q_points, + data.n_q_points_face); + } + + private: + template + static void + interpolate_generic(const Number * input, + Number * output, + const bool do_gradients, + const unsigned int face_no, + const unsigned int n_points_1d, + const AlignedVector *shape_data, + const unsigned int dofs_per_component_on_cell, + const unsigned int dofs_per_component_on_face) { internal::EvaluatorTensorProduct - evalf(data.data.front().shape_data_on_face[face_no % 2], + evalf(shape_data[face_no % 2], AlignedVector(), AlignedVector(), - data.data.front().fe_degree + 1, + n_points_1d, 0); - const unsigned int in_stride = do_evaluate ? - data.dofs_per_component_on_cell : - 2 * data.dofs_per_component_on_face; - const unsigned int out_stride = do_evaluate ? - 2 * data.dofs_per_component_on_face : - data.dofs_per_component_on_cell; + const unsigned int in_stride = + do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face; + const unsigned int out_stride = + do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell; const unsigned int face_direction = face_no / 2; for (unsigned int c = 0; c < n_components; c++) { diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index f20abbc98d..b30412de08 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -229,6 +229,18 @@ namespace internal */ AlignedVector shape_data_on_face[2]; + /** + * Collects all data of 1D nodal shape values (defined by the Lagrange + * polynomials in the points of the quadrature rule) evaluated at the + * point 0 and 1 (the vertices) in one data structure. + * + * This data structure can be used to interpolate from the cell to the + * face quadrature points. + * + * @note In contrast to shape_data_on_face, only the vales are evaluated. + */ + AlignedVector quadrature_data_on_face[2]; + /** * Stores one-dimensional values of shape functions on subface. Since * there are two subfaces, store two variants. diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 8b12453bab..41a90beb54 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -112,8 +112,10 @@ namespace internal univariate_shape_data.shape_gradients_collocation; auto &shape_hessians_collocation = univariate_shape_data.shape_hessians_collocation; - auto &inverse_shape_values = univariate_shape_data.inverse_shape_values; - auto &shape_data_on_face = univariate_shape_data.shape_data_on_face; + auto &inverse_shape_values = univariate_shape_data.inverse_shape_values; + auto &shape_data_on_face = univariate_shape_data.shape_data_on_face; + auto &quadrature_data_on_face = + univariate_shape_data.quadrature_data_on_face; auto &values_within_subface = univariate_shape_data.values_within_subface; auto &gradients_within_subface = univariate_shape_data.gradients_within_subface; @@ -298,6 +300,23 @@ namespace internal fe->shape_grad_grad(my_i, q_point)[0][0]; } + if (n_q_points_1d < 200) + { + quadrature_data_on_face[0].resize(quad.size() * 3); + quadrature_data_on_face[1].resize(quad.size() * 3); + + dealii::FE_DGQArbitraryNodes<1> fe_quad(quad); + + for (unsigned int i = 0; i < quad.size(); ++i) + { + Point<1> q_point; + q_point[0] = 0; + quadrature_data_on_face[0][i] = fe_quad.shape_value(i, q_point); + q_point[0] = 1; + quadrature_data_on_face[1][i] = fe_quad.shape_value(i, q_point); + } + } + // get gradient and Hessian transformation matrix for the polynomial // space associated with the quadrature rule (collocation space). We // need to avoid the case with more than a few hundreds of quadrature @@ -772,6 +791,8 @@ namespace internal { memory += MemoryConsumption::memory_consumption(shape_data_on_face[i]); + memory += + MemoryConsumption::memory_consumption(quadrature_data_on_face[i]); memory += MemoryConsumption::memory_consumption(values_within_subface[i]); memory += diff --git a/tests/matrix_free/interpolate_quadrature_01.cc b/tests/matrix_free/interpolate_quadrature_01.cc new file mode 100644 index 0000000000..8110ae3fa4 --- /dev/null +++ b/tests/matrix_free/interpolate_quadrature_01.cc @@ -0,0 +1,157 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + + +// Tests internal::FEFaceNormalEvaluationImpl::interpolate_quadrature() +// by comparing the results of FEFaceEvaluation::gather_evaluate(). + +template +class ExactSolution : public Function +{ +public: + ExactSolution() + {} + + virtual double + value(const Point &p, const unsigned int /*component*/ = 0) const + { + return p[0] + p[1]; + } +}; + +template > +void +test(const unsigned int n_refinements = 1) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria; + GridGenerator::hyper_cube(tria, 0, 1, true); + + tria.refine_global(n_refinements); + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ mapping(1); + QGauss<1> quad(n_points); + + AffineConstraints constraint; + + using MF = MatrixFree; + + typename MF::AdditionalData additional_data; + additional_data.mapping_update_flags = update_values; + additional_data.mapping_update_flags_inner_faces = update_values; + additional_data.mapping_update_flags_boundary_faces = update_values; + additional_data.mapping_update_flags_faces_by_cells = update_values; + additional_data.hold_all_faces_to_owned_cells = true; + + MF matrix_free; + matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); + + VectorType src, dst; + + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + VectorTools::interpolate(dof_handler, ExactSolution(), src); + + FEEvaluation phi( + matrix_free); + FEFaceEvaluation + phi_m(matrix_free, true); + + matrix_free.template loop_cell_centric( + [&](const auto &, auto &, const auto &src, const auto range) { + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + + phi.gather_evaluate(src, true, false); + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + face++) + { + phi_m.reinit(cell, face); + phi_m.gather_evaluate(src, true, false); + + AlignedVector temp( + phi_m.static_dofs_per_cell); + + internal::FEFaceNormalEvaluationImpl:: + template interpolate_quadrature( + matrix_free.get_shape_info(), + phi.begin_values(), + temp.data(), + false, + face); + + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + const auto u_cell = temp[q]; + const auto u_face = phi_m.get_value(q); + + for (unsigned int v = 0; v < VectorizedArray::size(); + ++v) + { + Assert(std::abs(u_cell[v] - u_face[v]) < 1e-10, + ExcMessage("Entries do not match!")); + } + } + } + } + }, + dst, + src); +} + +int +main() +{ + initlog(); + test<2, 1, 2, double, VectorizedArray>(); + test<3, 1, 2, double, VectorizedArray>(); + + deallog << "OK!" << std::endl; +} diff --git a/tests/matrix_free/interpolate_quadrature_01.output b/tests/matrix_free/interpolate_quadrature_01.output new file mode 100644 index 0000000000..5cfb783b8f --- /dev/null +++ b/tests/matrix_free/interpolate_quadrature_01.output @@ -0,0 +1,2 @@ + +DEAL::OK! -- 2.39.5