From ecde309bd8e493c7cb242efb0c45fb52bba1f352 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 13 Aug 2020 11:51:29 +0200 Subject: [PATCH] MF: Fix PBC in the case of non-standard orientation --- .../deal.II/matrix_free/face_setup_internal.h | 22 ++ tests/matrix_free/pbc_orientation_01.cc | 207 ++++++++++++++++++ tests/matrix_free/pbc_orientation_01.output | 9 + 3 files changed, 238 insertions(+) create mode 100644 tests/matrix_free/pbc_orientation_01.cc create mode 100644 tests/matrix_free/pbc_orientation_01.output diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h index 4f74facc14..bf154269a2 100644 --- a/include/deal.II/matrix_free/face_setup_internal.h +++ b/include/deal.II/matrix_free/face_setup_internal.h @@ -972,6 +972,28 @@ namespace internal cell->neighbor_of_coarser_neighbor(face_no).second; } + // special treatment of periodic boundaries + if (dim == 3 && cell->has_periodic_neighbor(face_no)) + { + const unsigned int exterior_face_orientation = + !cell->get_triangulation() + .get_periodic_face_map() + .at({cell, face_no}) + .second[0] + + 2 * cell->get_triangulation() + .get_periodic_face_map() + .at({cell, face_no}) + .second[1] + + 4 * cell->get_triangulation() + .get_periodic_face_map() + .at({cell, face_no}) + .second[2]; + + info.face_orientation = exterior_face_orientation; + + return info; + } + info.face_orientation = 0; const unsigned int interior_face_orientation = !cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) + diff --git a/tests/matrix_free/pbc_orientation_01.cc b/tests/matrix_free/pbc_orientation_01.cc new file mode 100644 index 0000000000..667eafbaf4 --- /dev/null +++ b/tests/matrix_free/pbc_orientation_01.cc @@ -0,0 +1,207 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// tests periodic boundaries for non-standard face orientations +// +// the grid is setup as in matrix_vector_faces_25 + +#include +#include + +#include + +#include + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + + +template +class ExactSolution : public Function +{ +public: + ExactSolution() + {} + + virtual double + value(const Point &p, const unsigned int /*component*/ = 0) const + { + // independent of p[2] since we apply PBC in z-direction + return p[0] + p[1]; + } +}; + + + +void generate_grid(Triangulation<3> &triangulation, int orientation) +{ + Point<3> vertices_1[] = {Point<3>(-0., -0., -0.), + Point<3>(+1., -0., -0.), + Point<3>(-0., +1., -0.), + Point<3>(+1., +1., -0.), + Point<3>(-0., -0., +0.5), + Point<3>(+1., -0., +0.5), + Point<3>(-0., +1., +0.5), + Point<3>(+1., +1., +0.5), + Point<3>(-0., -0., +1.), + Point<3>(+1., -0., +1.), + Point<3>(-0., +1., +1.), + Point<3>(+1., +1., +1.)}; + std::vector> vertices(&vertices_1[0], &vertices_1[12]); + + std::vector> cells(2, CellData<3>()); + + /* cell 0 */ + int cell_vertices_0[GeometryInfo<3>::vertices_per_cell] = { + 0, 1, 2, 3, 4, 5, 6, 7}; + + /* cell 1 */ + int cell_vertices_1[8][GeometryInfo<3>::vertices_per_cell] = { + {4, 5, 6, 7, 8, 9, 10, 11}, + {5, 7, 4, 6, 9, 11, 8, 10}, + {7, 6, 5, 4, 11, 10, 9, 8}, + {6, 4, 7, 5, 10, 8, 11, 9}, + {9, 8, 11, 10, 5, 4, 7, 6}, + {8, 10, 9, 11, 4, 6, 5, 7}, + {10, 11, 8, 9, 6, 7, 4, 5}, + {11, 9, 10, 8, 7, 5, 6, 4}}; + + for (const unsigned int j : GeometryInfo<3>::vertex_indices()) + { + cells[0].vertices[j] = cell_vertices_0[j]; + cells[1].vertices[j] = cell_vertices_1[orientation][j]; + } + cells[0].material_id = 0; + cells[1].material_id = 0; + + + triangulation.create_triangulation(vertices, cells, SubCellData()); +} + + + +template +void +test() +{ + if (dim == 2) + return; + + Triangulation<3> tria; + for (unsigned int orientation = 0; orientation < 8; ++orientation) + { + deallog << "Testing orientation case " << orientation << std::endl; + tria.clear(); + generate_grid(tria, orientation); + + FE_DGQ<3> fe(fe_degree); + DoFHandler<3> dof_handler(tria); + dof_handler.distribute_dofs(fe); + + for (auto &cell : tria) + for (auto &face : cell.face_iterators()) + { + if (std::abs(face->center()[2] - 0.0) < 10e-6) + face->set_boundary_id(4); + if (std::abs(face->center()[2] - 1.0) < 10e-6) + face->set_boundary_id(5); + } + + std::vector::cell_iterator>> + periodic_faces; + + dealii::GridTools::collect_periodic_faces(tria, 4, 5, 2, periodic_faces); + + tria.add_periodicity(periodic_faces); + + QGauss<1> quadrature(fe_degree + 1); + typename MatrixFree::AdditionalData additional_data; + additional_data.overlap_communication_computation = false; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points | + update_values); + additional_data.mapping_update_flags_inner_faces = + (update_JxW_values | update_normal_vectors | update_quadrature_points | + update_values); + additional_data.mapping_update_flags_boundary_faces = + (update_JxW_values | update_normal_vectors | update_quadrature_points | + update_values); + + MatrixFree data; + + AffineConstraints dummy; + dummy.close(); + data.reinit(dof_handler, dummy, quadrature, additional_data); + + using VectorType = LinearAlgebra::distributed::Vector; + + VectorType vec; + data.initialize_dof_vector(vec); + + VectorTools::interpolate(dof_handler, ExactSolution(), vec); + + + FEFaceEvaluation eval_minus( + data, true); + FEFaceEvaluation eval_plus( + data, false); + + data.template loop( + [](const auto &, auto &, const auto &, const auto) {}, + [&](const auto &, auto &, const auto &src, const auto face_range) { + for (unsigned int face = face_range.first; face < face_range.second; + face++) + { + eval_minus.reinit(face); + eval_minus.gather_evaluate(src, true, false); + eval_plus.reinit(face); + eval_plus.gather_evaluate(src, true, false); + + for (unsigned int q = 0; q < eval_minus.n_q_points; ++q) + { + const auto u_minus = eval_minus.get_value(q); + const auto u_plus = eval_plus.get_value(q); + + for (unsigned int v = 0; v < VectorizedArray::size(); + ++v) + { + Assert(std::abs(u_minus[v] - u_plus[v]) < 1e-10, + ExcMessage("Entries do not match!")); + } + } + } + }, + [](const auto &, auto &, const auto &, const auto) {}, + vec, + vec); + } +} + +int +main() +{ + initlog(); + test<3, 3>(); +} diff --git a/tests/matrix_free/pbc_orientation_01.output b/tests/matrix_free/pbc_orientation_01.output new file mode 100644 index 0000000000..4d7d8af862 --- /dev/null +++ b/tests/matrix_free/pbc_orientation_01.output @@ -0,0 +1,9 @@ + +DEAL::Testing orientation case 0 +DEAL::Testing orientation case 1 +DEAL::Testing orientation case 2 +DEAL::Testing orientation case 3 +DEAL::Testing orientation case 4 +DEAL::Testing orientation case 5 +DEAL::Testing orientation case 6 +DEAL::Testing orientation case 7 -- 2.39.5