From: Dominik Still Date: Mon, 15 Jul 2024 10:59:00 +0000 (+0200) Subject: Add face_to_cell_index_nodal for simplex elements X-Git-Tag: v9.6.0-rc1~90^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a682183e143414491b0c1df47a435a7e717754f8;p=dealii.git Add face_to_cell_index_nodal for simplex elements --- diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 9858856a05..8d12c435a7 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -496,6 +496,84 @@ namespace internal univariate_shape_data.nodal_at_cell_boundaries = true; + const ReferenceCell reference_cell = fe.reference_cell(); + if (reference_cell.is_simplex()) + { + if (dim == 2) + dofs_per_component_on_face = fe.degree + 1; + else + dofs_per_component_on_face = + (fe.degree + 1) * (fe.degree + 2) / 2; + + face_to_cell_index_nodal.reinit(reference_cell.n_faces(), + dofs_per_component_on_face); + + + for (unsigned int face = 0; face < reference_cell.n_faces(); + ++face) + { + // first get info from reference cell, i.e. the linear case + unsigned int d = 0; + for (; d < dim; ++d) + face_to_cell_index_nodal[face][d] = + reference_cell.face_to_cell_vertices(face, d, 1); + + // now fill the rest of the indices, start with the lines + if (fe.degree == 2) + for (; d < dofs_per_component_on_face; ++d) + face_to_cell_index_nodal[face][d] = + reference_cell.n_vertices() + + reference_cell.face_to_cell_lines(face, d - dim, 1); + + // in the cubic case it is more complicated as more DoFs are + // on the lines + else if (fe.degree == 3) + { + for (unsigned int line = 0; + d < dofs_per_component_on_face - 1; + ++line, d += 2) + { + const unsigned int face_to_cell_lines = + reference_cell.face_to_cell_lines(face, line, 1); + // check the direction of the line + // is it 0 -> 1 or 1 -> 0 + // as DoFs on the line are ordered differently + if (reference_cell.line_to_cell_vertices( + face_to_cell_lines, 0) == + reference_cell.face_to_cell_vertices(face, + line, + 1)) + { + face_to_cell_index_nodal[face][d] = + reference_cell.n_vertices() + + 2 * face_to_cell_lines; + face_to_cell_index_nodal[face][d + 1] = + reference_cell.n_vertices() + + 2 * face_to_cell_lines + 1; + } + else + { + face_to_cell_index_nodal[face][d + 1] = + reference_cell.n_vertices() + + 2 * face_to_cell_lines; + face_to_cell_index_nodal[face][d] = + reference_cell.n_vertices() + + 2 * face_to_cell_lines + 1; + } + } + // in 3D we also need the DoFs on the quads + if (dim == 3) + { + face_to_cell_index_nodal + [face][dofs_per_component_on_face - 1] = + reference_cell.n_vertices() + + 2 * reference_cell.n_lines() + face; + } + } + else if (fe.degree > 3) + DEAL_II_NOT_IMPLEMENTED(); + } + } // TODO: set up face_to_cell_index_nodal, face_to_cell_index_hermite, // face_orientations diff --git a/tests/matrix_free/face_to_cell_index_nodal_simplex.cc b/tests/matrix_free/face_to_cell_index_nodal_simplex.cc new file mode 100644 index 0000000000..1a8805dd0e --- /dev/null +++ b/tests/matrix_free/face_to_cell_index_nodal_simplex.cc @@ -0,0 +1,102 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2020 - 2024 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. +// +// ------------------------------------------------------------------------ + + + +// test the correctness of face_to_cell_index_nodal +// in internal::MatrixFreeFunctions::ShapeInfo +// for simplex elements + +#include + +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + + + +template +void +test(const FiniteElement &fe, const Quadrature &quad) +{ + for (unsigned int i = 0; i < fe.n_base_elements(); ++i) + { + internal::MatrixFreeFunctions::ShapeInfo shape_info; + shape_info.reinit(quad, fe, i); + deallog << "Testing: " << fe.get_name() << std::endl; + + Triangulation tria; + GridGenerator::reference_cell(tria, fe.reference_cell()); + + FE_SimplexP fe_continous(fe.degree); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe_continous); + + std::vector dof_indices( + fe_continous.n_dofs_per_face()); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + unsigned int face_counter = 0; + for (const auto &face : cell->face_iterators()) + { + face->get_dof_indices(dof_indices); + + unsigned int dof_counter = 0; + for (const auto i : dof_indices) + { + if (i == shape_info.face_to_cell_index_nodal[face_counter] + [dof_counter]) + deallog << "Ok "; + else + deallog << "Not ok"; + ++dof_counter; + } + deallog << std::endl; + ++face_counter; + } + } + } +} + + +int +main() +{ + initlog(); + for (unsigned int degree = 1; degree < 4; ++degree) + { + FE_SimplexP<2> fe(degree); + FE_SimplexDGP<2> fe_dg(degree); + QGaussSimplex<2> quad(1); + test(fe, quad); + test(fe_dg, quad); + + + FE_SimplexP<3> fe_3(degree); + FE_SimplexDGP<3> fe_dg_3(degree); + QGaussSimplex<3> quad_3(1); + test(fe_3, quad_3); + test(fe_dg_3, quad_3); + } +} diff --git a/tests/matrix_free/face_to_cell_index_nodal_simplex.output b/tests/matrix_free/face_to_cell_index_nodal_simplex.output new file mode 100644 index 0000000000..422b0e5d6b --- /dev/null +++ b/tests/matrix_free/face_to_cell_index_nodal_simplex.output @@ -0,0 +1,55 @@ + +DEAL::Testing: FE_SimplexP<2>(1) +DEAL::Ok Ok +DEAL::Ok Ok +DEAL::Ok Ok +DEAL::Testing: FE_SimplexDGP<2>(1) +DEAL::Ok Ok +DEAL::Ok Ok +DEAL::Ok Ok +DEAL::Testing: FE_SimplexP<3>(1) +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Testing: FE_SimplexDGP<3>(1) +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Testing: FE_SimplexP<2>(2) +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Testing: FE_SimplexDGP<2>(2) +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Ok Ok Ok +DEAL::Testing: FE_SimplexP<3>(2) +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Testing: FE_SimplexDGP<3>(2) +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok +DEAL::Testing: FE_SimplexP<2>(3) +DEAL::Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok +DEAL::Testing: FE_SimplexDGP<2>(3) +DEAL::Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok +DEAL::Testing: FE_SimplexP<3>(3) +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok +DEAL::Testing: FE_SimplexDGP<3>(3) +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok +DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok