From: Peter Munch Date: Thu, 31 Oct 2019 17:31:11 +0000 (+0100) Subject: Add new read_cell_data to FEEvaluationBase for arbitrary data types X-Git-Tag: v9.2.0-rc1~916^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F8981%2Fhead;p=dealii.git Add new read_cell_data to FEEvaluationBase for arbitrary data types --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 9716922d9d..0bb4baa57f 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -601,6 +601,16 @@ public: VectorizedArrayType read_cell_data(const AlignedVector &array) const; + /** + * The same as above, just for std::array of length of VectorizedArrayType for + * arbitrary data type. + */ + template + std::array + read_cell_data( + const AlignedVector> + &array) const; + //@} /** @@ -3578,6 +3588,40 @@ FEEvaluationBase:: +template +template +inline std::array +FEEvaluationBase:: + read_cell_data( + const AlignedVector> + &array) const +{ + Assert(matrix_info != nullptr, ExcNotImplemented()); + AssertDimension(array.size(), + matrix_info->get_task_info().cell_partition_data.back()); + if (is_face) + { + std::array out; + const unsigned int * cells = + is_interior_face ? + &this->matrix_info->get_face_info(cell).cells_interior[0] : + &this->matrix_info->get_face_info(cell).cells_exterior[0]; + for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i) + if (cells[i] != numbers::invalid_unsigned_int) + out[i] = array[cells[i] / VectorizedArrayType::n_array_elements] + [cells[i] % VectorizedArrayType::n_array_elements]; + return out; + } + else + return array[cell]; +} + + + namespace internal { // allows to select between block vectors and non-block vectors, which diff --git a/tests/matrix_free/read_cell_data.cc b/tests/matrix_free/read_cell_data.cc new file mode 100644 index 0000000000..415256267f --- /dev/null +++ b/tests/matrix_free/read_cell_data.cc @@ -0,0 +1,204 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 matrix-free read_cell_data for locally refined mesh + + +#include +#include + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template > +class Test +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + + void + run(unsigned int fe_degree) + { + Triangulation tria( + Triangulation::limit_level_difference_at_vertices); + GridGenerator::hyper_cube(tria); + + for (unsigned int i = 0; i < 3; i++) + { + tria.begin_active(i)->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + MappingQ mapping(1); + + QGauss<1> quad(fe_degree + 1); + + + + unsigned int level = numbers::invalid_unsigned_int; + { + AffineConstraints constraint; + + 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.mg_level = level; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraint, quad, additional_data); + + VectorType src, dst; + + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + this->setup_vector(matrix_free); + + matrix_free.loop(&Test::cell_operation, + &Test::face_operation, + &Test::boundary_operation, + this, + dst, + src); + + deallog << std::endl; + } + } + +private: + void + setup_vector(const MatrixFree &data) + { + unsigned int n_cells = data.n_cell_batches() + data.n_ghost_cell_batches(); + cell_ids.resize(n_cells); + + for (unsigned int cell = 0; cell < n_cells; cell++) + for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell); + lane++) + cell_ids[cell][lane] = data.get_cell_iterator(cell, lane)->id(); + } + + void + cell_operation(const MatrixFree &data, + VectorType &, + const VectorType &, + const std::pair &pair) const + { + FEEvaluation fe_eval(data); + for (auto cell = pair.first; cell < pair.second; cell++) + { + fe_eval.reinit(cell); + const auto cell_data = fe_eval.read_cell_data(cell_ids); + for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell); + lane++) + Assert(cell_data[lane] == data.get_cell_iterator(cell, lane)->id(), + ExcInternalError()); + } + } + + void + face_operation(const MatrixFree &data, + VectorType &, + const VectorType &, + const std::pair &pair) const + { + FEFaceEvaluation fe_eval_m(data, true); + FEFaceEvaluation fe_eval_p(data, false); + for (auto face = pair.first; face < pair.second; face++) + { + fe_eval_m.reinit(face); + fe_eval_p.reinit(face); + const auto cell_data_m = fe_eval_m.read_cell_data(cell_ids); + const auto cell_data_p = fe_eval_p.read_cell_data(cell_ids); + for (auto lane = 0u; lane < data.n_active_entries_per_face_batch(face); + lane++) + { + Assert(cell_data_m[lane] == + data.get_face_iterator(face, lane, true).first->id(), + ExcInternalError()); + Assert(cell_data_p[lane] == + data.get_face_iterator(face, lane, false).first->id(), + ExcInternalError()); + } + } + } + + + void + boundary_operation(const MatrixFree &data, + VectorType &, + const VectorType &, + const std::pair &pair) const + { + FEFaceEvaluation fe_eval(data, true); + for (auto face = pair.first; face < pair.second; face++) + { + fe_eval.reinit(face); + const auto cell_data = fe_eval.read_cell_data(cell_ids); + for (auto lane = 0u; lane < data.n_active_entries_per_face_batch(face); + lane++) + Assert(cell_data[lane] == + data.get_face_iterator(face, lane, true).first->id(), + ExcInternalError()); + } + } + + AlignedVector> + cell_ids; +}; + + +int +main() +{ + initlog(); + { + deallog.push("2d"); + Test<2> runner; + runner.run(1); + deallog.pop(); + } + { + deallog.push("3d"); + Test<3> runner; + runner.run(1); + deallog.pop(); + } +} diff --git a/tests/matrix_free/read_cell_data.output b/tests/matrix_free/read_cell_data.output new file mode 100644 index 0000000000..6ec9f053b6 --- /dev/null +++ b/tests/matrix_free/read_cell_data.output @@ -0,0 +1,3 @@ + +DEAL:2d:: +DEAL:3d::