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 <int dim,
+ typename Number,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+ 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<dim, Number, VectorizedArrayType> &matrix_free,
+ const AdditionalData &additional_data = AdditionalData())
+ {
+ this->matrix_free = &matrix_free;
+
+ std::vector<unsigned int> 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 <typename VectorTypeOut, typename VectorTypeIn>
+ void
+ cell_loop(const std::function<void(
+ const MatrixFree<dim, Number, VectorizedArrayType> &,
+ VectorTypeOut &,
+ const VectorTypeIn &,
+ const std::pair<unsigned int, unsigned int>)> &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<VectorTypeOut, VectorTypeIn>(
+ 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 <typename VectorTypeOut, typename VectorTypeIn>
+ void
+ loop(const std::function<
+ void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ VectorTypeOut &,
+ const VectorTypeIn &,
+ const std::pair<unsigned int, unsigned int>)> &cell_operation,
+ const std::function<
+ void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ VectorTypeOut &,
+ const VectorTypeIn &,
+ const std::pair<unsigned int, unsigned int>)> &face_operation,
+ const std::function<
+ void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ VectorTypeOut &,
+ const VectorTypeIn &,
+ const std::pair<unsigned int, unsigned int>,
+ 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<unsigned int>(category.first == fe_index_valid) +
+ static_cast<unsigned int>(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<VectorTypeOut, VectorTypeIn>(
+ 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<const MatrixFree<dim, Number, VectorizedArrayType>>
+ matrix_free;
+
+ /**
+ * Index of the valid FE. Currently, ony a single one is supported.
+ */
+ unsigned int fe_index_valid;
+ };
+
// implementations
#ifndef DOXYGEN
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+// Test MatrixFreeTools::ElementActivationAndDeactivationMatrixFree by comparing
+// the results with results obtained on fitted mesh.
+
+using namespace dealii;
+
+template <int dim>
+class Fu : public Function<dim>
+{
+public:
+ double
+ value(const Point<dim> &point, const unsigned int = 0) const
+ {
+ return std::sin(point[0] * 2 * numbers::PI) *
+ std::sin(point[1] * 2 * numbers::PI);
+ }
+};
+
+template <int dim>
+void
+test(const unsigned int n_refinements)
+{
+ using Number = double;
+ using VectorizedArrayType = VectorizedArray<Number>;
+ using VectorType = LinearAlgebra::distributed::Vector<Number>;
+ using FECellIntegrator =
+ FEEvaluation<dim, -1, 0, 1, Number, VectorizedArrayType>;
+ using FEFaceIntegrator =
+ FEFaceEvaluation<dim, -1, 0, 1, Number, VectorizedArrayType>;
+
+ for (unsigned int i = 0; i < 2; ++i)
+ {
+ parallel::distributed::Triangulation<dim> 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<dim> 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<dim>(FE_Q<dim>(1), FE_Nothing<dim>(1)));
+
+ typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+ data;
+ data.tasks_parallel_scheme =
+ MatrixFree<dim, double>::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<dim, Number, VectorizedArrayType> matrix_free;
+ matrix_free.reinit(MappingQ1<dim>(),
+ dof_handler,
+ AffineConstraints<Number>(),
+ QGauss<dim>(2),
+ data);
+
+ VectorType dst, src;
+ matrix_free.initialize_dof_vector(dst);
+ matrix_free.initialize_dof_vector(src);
+
+ VectorTools::interpolate(dof_handler, Fu<dim>(), 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<VectorType, VectorType>(
+ cell_operation, dst, src, true);
+
+ deallog << dst.l2_norm() << std::endl;
+
+ element_birth_and_death.template loop<VectorType, VectorType>(
+ cell_operation, face_operation, boundary_operation, dst, src, true);
+
+ deallog << dst.l2_norm() << std::endl;
+ }
+ else
+ {
+ matrix_free.template cell_loop<VectorType, VectorType>(cell_operation,
+ dst,
+ src,
+ true);
+
+ deallog << dst.l2_norm() << std::endl;
+
+ matrix_free.template loop<VectorType, VectorType>(
+ 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);
+}
projects / dealii.git / commitdiff