From 4156783ec270ee64b11108e9af88aba97e8e7488 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 15 Nov 2020 10:12:40 +0100 Subject: [PATCH] Make ShapeInfo compatible with Simplex::FE_P and Simplex::QGauss --- include/deal.II/matrix_free/shape_info.h | 7 +- .../matrix_free/shape_info.templates.h | 94 +++++++++ tests/simplex/matrix_free_shape_info_01.cc | 194 ++++++++++++++++++ ...ape_info_01.with_simplex_support=on.output | 1 + 4 files changed, 295 insertions(+), 1 deletion(-) create mode 100644 tests/simplex/matrix_free_shape_info_01.cc create mode 100644 tests/simplex/matrix_free_shape_info_01.with_simplex_support=on.output 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 }; diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 873300c9e5..802607bfd6 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -34,6 +34,8 @@ #include +#include + DEAL_II_NAMESPACE_OPEN @@ -134,7 +136,99 @@ namespace internal const FiniteElement &fe_in, const unsigned int base_element_number) { +#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT + if (quad_in.is_tensor_product() == false || + dynamic_cast *>(&fe_in) || + dynamic_cast *>(&fe_in)) + { + // specialization for arbitrary finite elements and quadrature rules + // as needed in the context, e.g., of simplices + + AssertDimension(dim, dim_q); + + const auto quad = Quadrature(quad_in); + const auto fe = &fe_in.base_element(base_element_number); + n_dimensions = dim; + n_components = fe_in.n_components(); + n_q_points = quad.size(); + dofs_per_component_on_cell = fe->n_dofs_per_cell(); + n_q_points_face = 0; // not implemented yet + dofs_per_component_on_face = 0; // + + Assert(fe->n_components() == 1, + ExcMessage( + "FEEvaluation only works for scalar finite elements.")); + + data.resize(1); + UnivariateShapeData &univariate_shape_data = data.front(); + data_access.reinit(n_dimensions, n_components); + data_access.fill(&univariate_shape_data); + + // note: we cannot write `univariate_shape_data.quadrature = quad`, + // since the quadrature rule within UnivariateShapeData expects + // a 1D quadrature rule. However, in this case we are not able to + // define that rule anyway so other code cannot use this information. + + univariate_shape_data.fe_degree = fe->degree; + univariate_shape_data.n_q_points_1d = quad.size(); + + if ((fe->n_dofs_per_cell() == 0) || (quad.size() == 0)) + return; + + // grant write access to common univariate shape data + auto &shape_values = univariate_shape_data.shape_values; + auto &shape_gradients = univariate_shape_data.shape_gradients; + + const unsigned int n_dofs = fe->n_dofs_per_cell(); + + const unsigned int array_size = n_dofs * n_q_points; + + shape_values.resize_fast(array_size); + shape_gradients.resize_fast(array_size * dim); + + for (unsigned int i = 0; i < n_dofs; ++i) + for (unsigned int q = 0; q < n_q_points; ++q) + { + shape_values[i * n_q_points + q] = + fe->shape_value(i, quad.point(q)); + + const auto grad = fe->shape_grad(i, quad.point(q)); + + for (int d = 0; d < dim; ++d) + shape_gradients[d * n_dofs * n_q_points + i * n_q_points + + q] = grad[d]; + } + + // TODO: also fill shape_hessians, inverse_shape_values, + // shape_data_on_face, quadrature_data_on_face, + // values_within_subface, gradients_within_subface, + // hessians_within_subface + + // note: shape_gradients_collocation, shape_hessians_collocation, + // shape_values_eo, shape_gradients_eo, shape_hessians_eo, + // shape_gradients_collocation_eo, shape_hessians_collocation_eo, + // inverse_shape_values_eo cannot be filled + + // indicate that no tensor product properties could be exploited + element_type = tensor_none; + univariate_shape_data.element_type = this->element_type; + + univariate_shape_data.nodal_at_cell_boundaries = true; + + lexicographic_numbering.resize(n_dofs); + + for (unsigned int i = 0; i < n_dofs; ++i) + lexicographic_numbering[i] = i; + + // TODO: setup face_to_cell_index_nodal, face_to_cell_index_hermite, + // face_orientations + + return; + } +#else Assert(quad_in.is_tensor_product(), ExcNotImplemented()); +#endif + const auto quad = quad_in.get_tensor_basis()[0]; const FiniteElement *fe = &fe_in.base_element(base_element_number); diff --git a/tests/simplex/matrix_free_shape_info_01.cc b/tests/simplex/matrix_free_shape_info_01.cc new file mode 100644 index 0000000000..cfe1180855 --- /dev/null +++ b/tests/simplex/matrix_free_shape_info_01.cc @@ -0,0 +1,194 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// Test ShapeData for Simplex::FE_P and Simplex::QGauss + +#include + +#include +#include + +#include + +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +class ExactSolution : public Function +{ +public: + ExactSolution() + {} + + virtual double + value(const Point &p, const unsigned int /*component*/ = 0) const + { + return p[0] + 2 * p[1]; + } +}; + +template +class FEEvaluationDummy +{ +public: + FEEvaluationDummy( + const internal::MatrixFreeFunctions::ShapeInfo &shape_info) + : shape_info(shape_info) + { + values.resize(shape_info.n_q_points); + gradients.resize(shape_info.n_q_points * dim); + } + + void + evaluate(const Vector &dof_values) + { + const unsigned int n_dofs = shape_info.dofs_per_component_on_cell; + const unsigned int n_q_points = shape_info.n_q_points; + + for (unsigned int q = 0, c = 0; q < shape_info.n_q_points; ++q) + { + values[q] = 0.0; + for (unsigned int i = 0; i < n_dofs; ++i, ++c) + values[q] += shape_info.data[0].shape_values[c] * dof_values[i]; + } + + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int q = 0; q < n_q_points; ++q) + { + gradients[n_q_points * d + q] = 0.0; + for (unsigned int i = 0; i < n_dofs; ++i) + gradients[n_q_points * d + q] += + shape_info.data[0] + .shape_gradients[d * n_dofs * n_q_points + i * n_q_points + q] * + dof_values[i]; + } + } + + Number + get_value(const unsigned int q_point) const + { + return values[q_point]; + } + + Tensor<1, dim, Number> + get_gradient(const unsigned int q_point) const + { + Tensor<1, dim, Number> out; + + for (unsigned int c = 0; c < dim; ++c) + out[c] = gradients[shape_info.n_q_points * c + q_point]; + + return out; + } + +private: + const internal::MatrixFreeFunctions::ShapeInfo &shape_info; + + AlignedVector values; + AlignedVector gradients; +}; + +template +void +test(const FiniteElement &fe) +{ + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingFE mapping(Simplex::FE_P(1)); + + Simplex::QGauss quadrature(1); + + internal::MatrixFreeFunctions::ShapeInfo shape_info(quadrature, fe); + + FEValues fe_values(mapping, + fe, + quadrature, + update_values | update_gradients | + update_inverse_jacobians); + + Vector src(dof_handler.n_dofs()); + VectorTools::interpolate(mapping, dof_handler, ExactSolution(), src); + + FEEvaluationDummy fe_eval(shape_info); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + + Vector dof_values(fe_values.dofs_per_cell); + cell->get_dof_values(src, dof_values); + fe_eval.evaluate(dof_values); + + // check values + { + std::vector values(fe_values.n_quadrature_points); + fe_values.get_function_values(src, values); + + for (const auto q : fe_values.quadrature_point_indices()) + { + Assert(std::abs(values[q] - fe_eval.get_value(q)) < 1e-10, + ExcMessage("Entries do not match!")); + } + } + + // check gradients + { + std::vector> gradients( + fe_values.n_quadrature_points); + fe_values.get_function_gradients(src, gradients); + + for (const auto q : fe_values.quadrature_point_indices()) + { + // apply the transposed inverse Jacobian of the mapping + Tensor<1, dim, Number> temp = fe_eval.get_gradient(q); + Tensor<1, dim, Number> temp_transformed; + for (int d = 0; d < dim; ++d) + { + Number sum = 0; + for (int e = 0; e < dim; ++e) + sum += fe_values.inverse_jacobian(q)[e][d] * temp[e]; + temp_transformed[d] = sum; + } + + for (int i = 0; i < dim; ++i) + Assert(std::abs(gradients[q][i] - temp_transformed[i]) < 1e-10, + ExcMessage("Entries do not match!")); + } + } + } +} + +int +main() +{ + initlog(); + + test<2>(Simplex::FE_P<2>(2)); +} diff --git a/tests/simplex/matrix_free_shape_info_01.with_simplex_support=on.output b/tests/simplex/matrix_free_shape_info_01.with_simplex_support=on.output new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/tests/simplex/matrix_free_shape_info_01.with_simplex_support=on.output @@ -0,0 +1 @@ + -- 2.39.5