From 897b217d7321b560f5367ebd869ecf52a80f632d Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Thu, 14 Jul 2022 10:20:10 +0200 Subject: [PATCH] Implement get_affine_constraints() --- include/deal.II/matrix_free/matrix_free.h | 17 ++++ .../matrix_free/matrix_free.templates.h | 45 ++++++++++ .../affine_constraints_store_and_get.cc | 84 +++++++++++++++++++ .../affine_constraints_store_and_get.output | 3 + 4 files changed, 149 insertions(+) create mode 100644 tests/matrix_free/affine_constraints_store_and_get.cc create mode 100644 tests/matrix_free/affine_constraints_store_and_get.output diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index d8e519ab6d..b98f7e32ac 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1630,6 +1630,15 @@ public: const DoFHandler & get_dof_handler(const unsigned int dof_handler_index = 0) const; + /** + * Return the AffineConstraints with the index as given to the + * respective `std::vector` argument in the reinit() function. Only available + * if the AffineConstraints objects have the same template parameter Number as + * MatrixFree. Throws an exception otherwise. + */ + const AffineConstraints & + get_affine_constraints(const unsigned int dof_handler_index = 0) const; + /** * Return the cell iterator in deal.II speak to a given cell batch * (populating several lanes in a VectorizedArray) and the lane index within @@ -2002,6 +2011,14 @@ private: */ std::vector>> dof_handlers; + /** + * Pointers to the AffineConstraints underlying the current problem. Only + * filled with an AffineConstraints object if objects of the same `Number` + * template parameter as the `Number` template of MatrixFree is passed to + * reinit(). Filled with nullptr otherwise. + */ + std::vector>> affine_constraints; + /** * Contains the information about degrees of freedom on the individual cells * and constraints. diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 90417b5f06..9919a0463e 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -228,6 +228,23 @@ MatrixFree::get_dof_handler( +template +const AffineConstraints & +MatrixFree::get_affine_constraints( + const unsigned int dof_handler_index) const +{ + AssertIndexRange(dof_handler_index, n_components()); + + AssertThrow(affine_constraints[dof_handler_index] != nullptr, + ExcMessage( + "AffineConstraints are only accessible if its Number type " + "matches with Number type of MatrixFree!")); + + return *(affine_constraints[dof_handler_index]); +} + + + template typename DoFHandler::cell_iterator MatrixFree::get_cell_iterator( @@ -330,6 +347,29 @@ MatrixFree::copy_from( +namespace internal +{ + template + void + store_affine_constraints( + const dealii::AffineConstraints *, + SmartPointer> &stored_constraints) + { + stored_constraints = nullptr; + } + + template + void + store_affine_constraints( + const dealii::AffineConstraints * affine_constraints, + SmartPointer> &stored_constraints) + { + stored_constraints = affine_constraints; + } +} // namespace internal + + + template template void @@ -374,6 +414,11 @@ MatrixFree::internal_reinit( .reinit(quad[nq][q_no], dof_handler[no]->get_fe(fe_no), b); } + // Store pointers to AffineConstraints objects if Number type matches + affine_constraints.resize(constraints.size()); + for (unsigned int no = 0; no < constraints.size(); ++no) + internal::store_affine_constraints(constraints[no], affine_constraints[no]); + if (additional_data.initialize_indices == true) { clear(); diff --git a/tests/matrix_free/affine_constraints_store_and_get.cc b/tests/matrix_free/affine_constraints_store_and_get.cc new file mode 100644 index 0000000000..c0d96d58b9 --- /dev/null +++ b/tests/matrix_free/affine_constraints_store_and_get.cc @@ -0,0 +1,84 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 - 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. +// +// --------------------------------------------------------------------- + + + +// test internal::store_affine_constraints() and get_affine_constraints() + +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + + +int +main() +{ + initlog(); + + const unsigned int fe_degree = 1; + constexpr int dim = 2; + + Triangulation tria; + GridGenerator::hyper_cube(tria, 0, 1, true); + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + QGauss<1> quadrature(fe_degree + 1); + + AffineConstraints double_ac; + double_ac.close(); + + { + typename MatrixFree::AdditionalData additional_data_double; + MatrixFree matrix_free_double; + matrix_free_double.reinit(MappingQ1{}, + dof_handler, + double_ac, + quadrature, + additional_data_double); + + if (&matrix_free_double.get_affine_constraints() == &double_ac) + deallog << "OK" << std::endl; + } + + { + typename MatrixFree::AdditionalData additional_data_float; + MatrixFree matrix_free_float; + matrix_free_float.reinit(MappingQ1{}, + dof_handler, + double_ac, + quadrature, + additional_data_float); + + try + { + matrix_free_float.get_affine_constraints(); + } + catch (...) + { + deallog << "OK" << std::endl; + } + } +} diff --git a/tests/matrix_free/affine_constraints_store_and_get.output b/tests/matrix_free/affine_constraints_store_and_get.output new file mode 100644 index 0000000000..8b3b075900 --- /dev/null +++ b/tests/matrix_free/affine_constraints_store_and_get.output @@ -0,0 +1,3 @@ + +DEAL::OK +DEAL::OK -- 2.39.5