From: Peter Munch Date: Tue, 1 Sep 2020 10:50:06 +0000 (+0200) Subject: Add MatrixFreeTools::categorize_accoring_boundary_ids_for_ecl X-Git-Tag: v9.3.0-rc1~1132^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10871%2Fhead;p=dealii.git Add MatrixFreeTools::categorize_accoring_boundary_ids_for_ecl --- diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index cfffa41214..ffcdd5b790 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -18,6 +18,7 @@ #include +#include DEAL_II_NAMESPACE_OPEN @@ -27,7 +28,83 @@ DEAL_II_NAMESPACE_OPEN * evaluation. */ namespace MatrixFreeTools -{} // namespace MatrixFreeTools +{ + /** + * Modify @p additional_data so that cells are categorized + * according to their boundary IDs, making face integrals in the case of + * cell-centric loop simpler. + */ + template + void + categorize_by_boundary_ids( + const Triangulation &tria, + AdditionalData & additional_data, + const unsigned int level = numbers::invalid_unsigned_int) + { + bool is_mg = (level != numbers::invalid_unsigned_int); + + // ... create list for the category of each cell + if (is_mg) + additional_data.cell_vectorization_category.resize( + std::distance(tria.begin(level), tria.end(level))); + else + additional_data.cell_vectorization_category.resize(tria.n_active_cells()); + + // ... set up scaling factor + std::vector factors(dim * 2); + + const auto bids = tria.get_boundary_ids(); + + std::map bid_map; + for (const auto id : tria.get_boundary_ids()) + bid_map[id] = bid_map.size() + 1; + + { + unsigned int n_bids = bids.size() + 1; + int offset = 1; + for (unsigned int i = 0; i < dim * 2; i++, offset = offset * n_bids) + factors[i] = offset; + } + + auto to_category = [&](auto &cell) { + unsigned int c_num = 0; + for (unsigned int i = 0; i < dim * 2; i++) + { + auto &face = *cell->face(i); + if (face.at_boundary()) + c_num += factors[i] * bid_map[face.boundary_id()]; + } + return c_num; + }; + + if (!is_mg) + { + for (auto cell = tria.begin_active(); cell != tria.end(); ++cell) + { + if (cell->is_locally_owned()) + additional_data + .cell_vectorization_category[cell->active_cell_index()] = + to_category(cell); + } + } + else + { + for (auto cell = tria.begin(level); cell != tria.end(level); ++cell) + { + if (cell->is_locally_owned_on_level()) + additional_data.cell_vectorization_category[cell->index()] = + to_category(cell); + } + } + + // ... finalize set up of matrix_free + additional_data.hold_all_faces_to_owned_cells = true; + additional_data.cell_vectorization_categories_strict = true; + additional_data.mapping_update_flags_faces_by_cells = + additional_data.mapping_update_flags_inner_faces | + additional_data.mapping_update_flags_boundary_faces; + } +} // namespace MatrixFreeTools DEAL_II_NAMESPACE_CLOSE diff --git a/tests/matrix_free/categorize_by_boundary_ids_01.cc b/tests/matrix_free/categorize_by_boundary_ids_01.cc new file mode 100644 index 0000000000..e0e582a45b --- /dev/null +++ b/tests/matrix_free/categorize_by_boundary_ids_01.cc @@ -0,0 +1,128 @@ +// --------------------------------------------------------------------- +// +// 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 MatrixFreeTools::categorize_by_boundary_ids() + +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + + +template +void +test(const unsigned int n_refinements) +{ + const unsigned int fe_degree = 1; + + Triangulation tria; + GridGenerator::hyper_cube(tria, 0, 1, true); + + tria.refine_global(n_refinements); + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + 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.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + AffineConstraints dummy; + dummy.close(); + + const auto process = [&]() { + bool result = true; + + MatrixFree data; + data.reinit(dof_handler, dummy, quadrature, additional_data); + + using VectorType = Vector; + + VectorType vec; + data.initialize_dof_vector(vec); + + data.template cell_loop( + [&](const auto &, auto &, const auto &, const auto cell_range) { + for (unsigned int cell = cell_range.first; cell < cell_range.second; + ++cell) + { + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + bool temp = true; + + for (unsigned int v = 1; v < data.n_components_filled(cell); + ++v) + temp &= (data.get_faces_by_cells_boundary_id(cell, face)[0] == + data.get_faces_by_cells_boundary_id(cell, face)[v]); + + result &= temp; + } + } + }, + vec, + vec); + + return result; + }; + + AssertDimension(process(), false); // not categorize + + MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data); + + AssertDimension(process(), true); // categorized + + deallog << "OK!" << std::endl; +} + +int +main() +{ + initlog(); + + { + deallog.push("2D"); + test<2>(3); + deallog.pop(); + } + + { + deallog.push("3D"); + test<3>(3); + deallog.pop(); + } +} diff --git a/tests/matrix_free/categorize_by_boundary_ids_01.output b/tests/matrix_free/categorize_by_boundary_ids_01.output new file mode 100644 index 0000000000..417aae8b90 --- /dev/null +++ b/tests/matrix_free/categorize_by_boundary_ids_01.output @@ -0,0 +1,3 @@ + +DEAL:2D::OK! +DEAL:3D::OK!