From e931728ecc42b6162c6afeda52a6795bf445c5a9 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 27 Nov 2020 09:08:58 +0100 Subject: [PATCH] MatrixFree: set up shape_values_face and shape_gradients_face --- .../matrix_free/mapping_info.templates.h | 25 ++++--- .../matrix_free/matrix_free.templates.h | 15 ++-- include/deal.II/matrix_free/shape_info.h | 14 ++++ .../matrix_free/shape_info.templates.h | 69 ++++++++++++++++++- include/deal.II/matrix_free/util.h | 53 ++++++++++++++ 5 files changed, 160 insertions(+), 16 deletions(-) create mode 100644 include/deal.II/matrix_free/util.h diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 984795d3c2..444d5c1b8c 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -30,6 +30,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -71,7 +72,7 @@ namespace internal for (unsigned int i = 0; i < n_q_points; ++i) quadrature_weights[i] = quadrature.weight(i); - // note: quadrature_1 and tensor_quadrature_weights are not set up + // note: quadrature_1d and tensor_quadrature_weights are not set up // TODO: set up face_orientations (void)update_flags_inner_faces; @@ -404,15 +405,17 @@ namespace internal if (flag == false) { #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT - continue; // TODO: face data is not set up + const auto quad_face = get_face_quadrature(quad[my_q][hpq]); + face_data[my_q].descriptor[hpq].initialize(quad_face, + update_default); #else Assert(false, ExcNotImplemented()); #endif } - - face_data[my_q].descriptor[hpq].initialize( - quad[my_q][hpq].get_tensor_basis()[0], - update_flags_boundary_faces); + else + face_data[my_q].descriptor[hpq].initialize( + quad[my_q][hpq].get_tensor_basis()[0], + update_flags_boundary_faces); } face_data_by_cells[my_q].descriptor.resize(n_hp_quads); @@ -428,14 +431,16 @@ namespace internal if (flag == false) { #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT - continue; // TODO: face data is not set up + const auto quad_face = get_face_quadrature(quad[my_q][hpq]); + face_data_by_cells[my_q].descriptor[hpq].initialize( + quad_face, update_default); #else Assert(false, ExcNotImplemented()); #endif } - - face_data_by_cells[my_q].descriptor[hpq].initialize( - quad[my_q][hpq].get_tensor_basis()[0], update_default); + else + face_data_by_cells[my_q].descriptor[hpq].initialize( + quad[my_q][hpq].get_tensor_basis()[0], update_default); } } diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index f2e2325cf1..e2c0e006f3 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -41,6 +41,8 @@ #include #include +#include + #ifdef DEAL_II_WITH_TBB # include @@ -1335,7 +1337,8 @@ MatrixFree::initialize_indices( Table<2, internal::MatrixFreeFunctions::ShapeInfo> shape_info_dummy( shape_info.size(0), shape_info.size(2)); { - QGauss quad(1); + Quadrature quad(QGauss(1)); + Quadrature quad_simplex(Simplex::QGauss(1)); for (unsigned int no = 0, c = 0; no < dof_handlers.size(); no++) for (unsigned int b = 0; b < dof_handlers[no]->get_fe(0).n_base_elements(); @@ -1343,9 +1346,13 @@ MatrixFree::initialize_indices( for (unsigned int fe_no = 0; fe_no < dof_handlers[no]->get_fe_collection().size(); ++fe_no) - shape_info_dummy(c, fe_no).reinit(quad, - dof_handlers[no]->get_fe(fe_no), - b); + shape_info_dummy(c, fe_no).reinit( + dof_handlers[no]->get_fe(fe_no).reference_cell_type() == + ReferenceCell::get_hypercube(dim) ? + quad : + quad_simplex, + dof_handlers[no]->get_fe(fe_no), + b); } const unsigned int n_lanes = VectorizedArrayType::size(); diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 8b7eba41b8..23d07de2bb 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -282,6 +282,20 @@ namespace internal * end points of the unit cell. */ bool nodal_at_cell_boundaries; + + /** + * Stores the shape values of the finite element evaluated on all + * quadrature points for all faces and orientations (no tensor-product + * structure exploited). + */ + Table<3, Number> shape_values_face; + + /** + * Stores the shape gradients of the finite element evaluated on all + * quadrature points for all faces, orientations, and directions + * (no tensor-product structure exploited). + */ + Table<4, Number> shape_gradients_face; }; diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 81280e4f48..90f7cea9c3 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -30,9 +31,12 @@ #include #include +#include + #include #include +#include #include @@ -184,8 +188,11 @@ namespace internal 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; + auto &shape_values = univariate_shape_data.shape_values; + auto &shape_values_face = univariate_shape_data.shape_values_face; + auto &shape_gradients = univariate_shape_data.shape_gradients; + auto &shape_gradients_face = + univariate_shape_data.shape_gradients_face; const unsigned int n_dofs = fe->n_dofs_per_cell(); @@ -207,6 +214,64 @@ namespace internal q] = grad[d]; } + const auto reference_cell_type = ReferenceCell::get_simplex(dim); + + try + { + const auto quad_face = get_face_quadrature(quad); + this->n_q_points_face = quad_face.size(); + + const unsigned int n_face_orientations = dim == 2 ? 2 : 6; + const unsigned int n_faces = + ReferenceCell::internal::Info::get_cell(reference_cell_type) + .n_faces(); + + const auto projected_quad_face = + QProjector::project_to_all_faces(reference_cell_type, + quad_face); + + shape_values_face.reinit( + {n_faces, n_face_orientations, n_dofs * n_q_points_face}); + + shape_gradients_face.reinit( + {n_faces, n_face_orientations, dim, n_dofs * n_q_points_face}); + + for (unsigned int f = 0; f < n_faces; ++f) + for (unsigned int o = 0; o < n_face_orientations; ++o) + { + const auto offset = + QProjector::DataSetDescriptor::face( + reference_cell_type, + f, + (o ^ 1) & 1, // face_orientation + (o >> 1) & 1, // face_flip + (o >> 2) & 1, // face_rotation + n_q_points_face); + + for (unsigned int i = 0; i < n_dofs; ++i) + for (unsigned int q = 0; q < n_q_points_face; ++q) + { + const auto point = + projected_quad_face.point(q + offset); + + shape_values_face(f, o, i * n_q_points_face + q) = + fe->shape_value(i, point); + + const auto grad = fe->shape_grad(i, point); + + for (int d = 0; d < dim; ++d) + shape_gradients_face( + f, o, d, i * n_q_points_face + q) = grad[d]; + } + } + } + catch (...) + { + // TODO: this might happen if the quadrature rule and the + // the FE do not match + this->n_q_points_face = 0; + } + // TODO: also fill shape_hessians, inverse_shape_values, // shape_data_on_face, quadrature_data_on_face, // values_within_subface, gradients_within_subface, diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h new file mode 100644 index 0000000000..b46f3c27f5 --- /dev/null +++ b/include/deal.II/matrix_free/util.h @@ -0,0 +1,53 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_util_h +#define dealii_matrix_free_util_h + + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + namespace MatrixFreeFunctions + { + template + inline Quadrature + get_face_quadrature(const Quadrature &quad) + { + if (dim == 2 || dim == 3) + for (unsigned int i = 1; i <= 3; ++i) + if (quad == Simplex::QGauss(i)) + return Simplex::QGauss(i); + + AssertThrow(false, ExcNotImplemented()); + + return Quadrature(); + } + + } // end of namespace MatrixFreeFunctions +} // end of namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif -- 2.39.5