From: Peter Munch Date: Sat, 9 Jul 2022 08:37:56 +0000 (+0200) Subject: Add MatrixFreeTools::ElementActivationAndDeactivationMatrixFree X-Git-Tag: v9.5.0-rc1~994^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14121%2Fhead;p=dealii.git Add MatrixFreeTools::ElementActivationAndDeactivationMatrixFree --- diff --git a/doc/news/changes/minor/20220709Munch b/doc/news/changes/minor/20220709Munch new file mode 100644 index 0000000000..3160896c0a --- /dev/null +++ b/doc/news/changes/minor/20220709Munch @@ -0,0 +1,7 @@ +New: The new class MatrixFreeTools::ElementActivationAndDeactivationMatrixFree +is a wrapper around MatrixFree designed to deal with +DoFHandler objects involving cells without degrees of freedom, i.e., +cells using FENothing as element type. The class helps to implement the +"element birth and death technique". +
+(Peter Munch, Sebastian Proell, 2022/07/09) diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index c5b34249d0..61082de1bd 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -157,6 +157,188 @@ namespace MatrixFreeTools const unsigned int first_selected_component = 0); + + /** + * A wrapper around MatrixFree to help users to deal with DoFHandler + * objects involving cells without degrees of freedom, i.e., + * cells using FENothing as element type. In the following we call such cells + * deactivated. All other cells are activated. In contrast to MatrixFree, + * this class skips deactivated cells and faces between activated and + * deactivated cells are treated as boundary faces. + */ + template > + class ElementActivationAndDeactivationMatrixFree + { + public: + /** + * Struct that helps to configure + * ElementActivationAndDeactivationMatrixFree. + */ + struct AdditionalData + { + /** + * Constructor. + */ + AdditionalData(const unsigned int dof_index = 0) + : dof_index(dof_index) + {} + + /** + * Index of the DoFHandler within MatrixFree to be used. + */ + unsigned int dof_index; + }; + + /** + * Reinitialize class based on a given MatrixFree instance. + * Particularly, the index of the valid FE is determined. + * + * @note At the moment, only DoFHandler objects are accepted + * that are initialized with FECollection objects with at most + * two finite elements (i.e., `FENothing` and another finite + * element). + */ + void + reinit(const MatrixFree &matrix_free, + const AdditionalData &additional_data = AdditionalData()) + { + this->matrix_free = &matrix_free; + + std::vector valid_fe_indices; + + const auto &fe_collection = + matrix_free.get_dof_handler(additional_data.dof_index) + .get_fe_collection(); + + for (unsigned int i = 0; i < fe_collection.size(); ++i) + if (fe_collection[i].n_dofs_per_cell() > 0) + valid_fe_indices.push_back(i); + + // TODO: relax this so that arbitrary number of valid + // FEs can be accepted + AssertDimension(valid_fe_indices.size(), 1); + + fe_index_valid = *valid_fe_indices.begin(); + } + + /** + * Loop over all activated cells. + * + * For the meaning of the parameters see MatrixFree::cell_loop(). + */ + template + void + cell_loop(const std::function &, + VectorTypeOut &, + const VectorTypeIn &, + const std::pair)> &cell_operation, + VectorTypeOut & dst, + const VectorTypeIn & src, + const bool zero_dst_vector = false) + { + const auto ebd_cell_operation = [&](const auto &matrix_free, + auto & dst, + const auto &src, + const auto range) { + const auto category = matrix_free.get_cell_range_category(range); + + if (category != fe_index_valid) + return; + + cell_operation(matrix_free, dst, src, range); + }; + + matrix_free->template cell_loop( + ebd_cell_operation, dst, src, zero_dst_vector); + } + + /** + * Loop over all activated cells and faces. Faces between activated and + * deactivated cells are treated as boundary faces with the boundary ID + * numbers::internal_face_boundary_id. + * + * For the meaning of the parameters see MatrixFree::cell_loop(). + */ + template + void + loop(const std::function< + void(const MatrixFree &, + VectorTypeOut &, + const VectorTypeIn &, + const std::pair)> &cell_operation, + const std::function< + void(const MatrixFree &, + VectorTypeOut &, + const VectorTypeIn &, + const std::pair)> &face_operation, + const std::function< + void(const MatrixFree &, + VectorTypeOut &, + const VectorTypeIn &, + const std::pair, + const bool)> &boundary_operation, + VectorTypeOut & dst, + const VectorTypeIn & src, + const bool zero_dst_vector = false) + { + const auto ebd_cell_operation = [&](const auto &matrix_free, + auto & dst, + const auto &src, + const auto range) { + const auto category = matrix_free.get_cell_range_category(range); + + if (category != fe_index_valid) + return; + + cell_operation(matrix_free, dst, src, range); + }; + + const auto ebd_internal_or_boundary_face_operation = + [&](const auto &matrix_free, + auto & dst, + const auto &src, + const auto range) { + const auto category = matrix_free.get_face_range_category(range); + + const unsigned int type = + static_cast(category.first == fe_index_valid) + + static_cast(category.second == fe_index_valid); + + if (type == 0) + return; // deactivated face -> nothing to do + + if (type == 1) // boundary face + boundary_operation( + matrix_free, dst, src, range, category.first == fe_index_valid); + else if (type == 2) // internal face + face_operation(matrix_free, dst, src, range); + }; + + matrix_free->template loop( + ebd_cell_operation, + ebd_internal_or_boundary_face_operation, + ebd_internal_or_boundary_face_operation, + dst, + src, + zero_dst_vector); + } + + private: + /** + * Reference to the underlying MatrixFree object. + */ + SmartPointer> + matrix_free; + + /** + * Index of the valid FE. Currently, ony a single one is supported. + */ + unsigned int fe_index_valid; + }; + // implementations #ifndef DOXYGEN diff --git a/tests/matrix_free/element_birth_and_death_01.cc b/tests/matrix_free/element_birth_and_death_01.cc new file mode 100644 index 0000000000..53171d0366 --- /dev/null +++ b/tests/matrix_free/element_birth_and_death_01.cc @@ -0,0 +1,209 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +// Test MatrixFreeTools::ElementActivationAndDeactivationMatrixFree by comparing +// the results with results obtained on fitted mesh. + +using namespace dealii; + +template +class Fu : public Function +{ +public: + double + value(const Point &point, const unsigned int = 0) const + { + return std::sin(point[0] * 2 * numbers::PI) * + std::sin(point[1] * 2 * numbers::PI); + } +}; + +template +void +test(const unsigned int n_refinements) +{ + using Number = double; + using VectorizedArrayType = VectorizedArray; + using VectorType = LinearAlgebra::distributed::Vector; + using FECellIntegrator = + FEEvaluation; + using FEFaceIntegrator = + FEFaceEvaluation; + + for (unsigned int i = 0; i < 2; ++i) + { + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + if (i == 0) + GridGenerator::subdivided_hyper_rectangle( + tria, {2, 2}, {0.0, 0.0}, {1.0, 1.0}, true); + else + GridGenerator::subdivided_hyper_rectangle( + tria, {2, 1}, {0.0, 0.0}, {1.0, 0.5}, true); + + tria.refine_global(n_refinements); + + DoFHandler dof_handler(tria); + + for (const auto &cell : + filter_iterators(dof_handler.active_cell_iterators(), + IteratorFilters::LocallyOwnedCell())) + { + if (cell->center()[1] < 0.5) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); + } + + dof_handler.distribute_dofs( + hp::FECollection(FE_Q(1), FE_Nothing(1))); + + typename MatrixFree::AdditionalData + data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients; + data.mapping_update_flags_boundary_faces = update_gradients; + data.mapping_update_flags_inner_faces = update_gradients; + + MatrixFree matrix_free; + matrix_free.reinit(MappingQ1(), + dof_handler, + AffineConstraints(), + QGauss(2), + data); + + VectorType dst, src; + matrix_free.initialize_dof_vector(dst); + matrix_free.initialize_dof_vector(src); + + VectorTools::interpolate(dof_handler, Fu(), src); + + const auto cell_operation = [&](const auto &matrix_free, + auto & dst, + const auto &src, + const auto range) { + FECellIntegrator phi(matrix_free); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(src, EvaluationFlags::gradients); + + for (const auto q : phi.quadrature_point_indices()) + phi.submit_gradient(phi.get_gradient(q), q); + + phi.integrate_scatter(EvaluationFlags::gradients, dst); + } + }; + + const auto face_operation = + [&](const auto &, auto &, const auto &, const auto) { + // nothing to do here but in the DG case + }; + + const auto boundary_operation = [&](const auto &matrix_free, + auto &, + const auto &, + const auto range, + const bool is_interior_face = true) { + FEFaceIntegrator phi(matrix_free, is_interior_face); + + for (unsigned int face = range.first; face < range.second; ++face) + { + phi.reinit(face); + + const double scaling = + phi.at_boundary() ? (phi.boundary_id() + 1.0) : (2.0 * dim); + + phi.gather_evaluate(src, EvaluationFlags::gradients); + + for (const auto q : phi.quadrature_point_indices()) + { + auto gradient = phi.get_gradient(q); + auto normal = phi.get_normal_vector(q); + + if (is_interior_face == false) // fix sign! + normal *= -1.0; + + phi.submit_value(gradient * normal * scaling, q); + } + + phi.integrate_scatter(EvaluationFlags::values, dst); + } + }; + + if (i == 0) + { + MatrixFreeTools::ElementActivationAndDeactivationMatrixFree< + dim, + Number, + VectorizedArrayType> + element_birth_and_death; + + element_birth_and_death.reinit(matrix_free); + + element_birth_and_death.template cell_loop( + cell_operation, dst, src, true); + + deallog << dst.l2_norm() << std::endl; + + element_birth_and_death.template loop( + cell_operation, face_operation, boundary_operation, dst, src, true); + + deallog << dst.l2_norm() << std::endl; + } + else + { + matrix_free.template cell_loop(cell_operation, + dst, + src, + true); + + deallog << dst.l2_norm() << std::endl; + + matrix_free.template loop( + cell_operation, face_operation, boundary_operation, dst, src, true); + + deallog << dst.l2_norm() << std::endl; + } + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + mpi_initlog(); + test<2>(2); +} diff --git a/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=1.output b/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..b83d53d551 --- /dev/null +++ b/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=1.output @@ -0,0 +1,5 @@ + +DEAL::3.59694 +DEAL::9.14748 +DEAL::3.59694 +DEAL::9.14748 diff --git a/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=2.output b/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..b83d53d551 --- /dev/null +++ b/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=2.output @@ -0,0 +1,5 @@ + +DEAL::3.59694 +DEAL::9.14748 +DEAL::3.59694 +DEAL::9.14748