From 85e1360c034802bce0be24641517a50abb55fe90 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 10 Dec 2021 16:46:44 +0100 Subject: [PATCH] MatrixFree: fix orientation for cell-centric loop and pbc --- .../deal.II/matrix_free/evaluation_kernels.h | 2 + include/deal.II/matrix_free/fe_evaluation.h | 2 +- tests/matrix_free/pbc_orientation_02.cc | 217 ++++++++++++++++++ tests/matrix_free/pbc_orientation_02.output | 17 ++ 4 files changed, 237 insertions(+), 1 deletion(-) create mode 100644 tests/matrix_free/pbc_orientation_02.cc create mode 100644 tests/matrix_free/pbc_orientation_02.output diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index dc84931d63..81ff490f76 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3996,6 +3996,8 @@ namespace internal proc.value( temp1[reorientate(v, i)][v], vector_ptrs[v][index_array_nodal[v][i]]); + else if (integrate == false) + temp1[i][v] = 0.0; } } } diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index be9a88ce8f..b171e4344b 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -8919,7 +8919,7 @@ FEFaceEvaluation table{ - {0, 1, 0, 3, 6, 5, 4, 7}}; + {0, 1, 2, 3, 6, 5, 4, 7}}; face_orientation = table[face_orientation]; } diff --git a/tests/matrix_free/pbc_orientation_02.cc b/tests/matrix_free/pbc_orientation_02.cc new file mode 100644 index 0000000000..1c1b63518f --- /dev/null +++ b/tests/matrix_free/pbc_orientation_02.cc @@ -0,0 +1,217 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// Like pbc_orientation_02 but with cell-centric loops. + +#include +#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[orientation < 8 ? 0 : 1].vertices[j] = cell_vertices_0[j]; + cells[orientation < 8 ? 1 : 0].vertices[j] = + cell_vertices_1[orientation % 8][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 < 16; ++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.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.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + MatrixFreeTools::categorize_by_boundary_ids( + dof_handler.get_triangulation(), additional_data); + + 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); + + + + data.template loop_cell_centric( + [](const auto &data, auto &, const auto &src, const auto cell_range) { + FEFaceEvaluation eval_minus( + data, true); + FEFaceEvaluation eval_plus( + data, false); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; + ++cell) + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + const auto boundary_ids = + data.get_faces_by_cells_boundary_id(cell, face); + const auto boundary_id = boundary_ids[0]; + + if (boundary_id != numbers::internal_face_boundary_id) + continue; + + eval_minus.reinit(cell, face); + eval_minus.gather_evaluate(src, true, false); + eval_plus.reinit(cell, 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!")); + } + } + } + }, + vec, + vec); + } +} + +int +main() +{ + initlog(); + test<3, 3>(); +} diff --git a/tests/matrix_free/pbc_orientation_02.output b/tests/matrix_free/pbc_orientation_02.output new file mode 100644 index 0000000000..c822334ac8 --- /dev/null +++ b/tests/matrix_free/pbc_orientation_02.output @@ -0,0 +1,17 @@ + +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 +DEAL::Testing orientation case 8 +DEAL::Testing orientation case 9 +DEAL::Testing orientation case 10 +DEAL::Testing orientation case 11 +DEAL::Testing orientation case 12 +DEAL::Testing orientation case 13 +DEAL::Testing orientation case 14 +DEAL::Testing orientation case 15 -- 2.39.5