From f29e40400d6a7a82417f3a9e5f8846ec36bdc1b9 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 17 Feb 2023 15:35:29 -0500 Subject: [PATCH] implement VectorTools::add_constant --- doc/news/changes/minor/20230217tjhei | 4 + .../numerics/vector_tools_mean_value.h | 32 +++ .../vector_tools_mean_value.templates.h | 149 +++++++++++++ .../numerics/vector_tools_mean_value.inst.in | 6 + tests/vector_tools/add_constant_01.cc | 176 +++++++++++++++ tests/vector_tools/add_constant_01.output | 2 + tests/vector_tools/add_constant_02.cc | 203 ++++++++++++++++++ .../add_constant_02.mpirun=2.output | 5 + 8 files changed, 577 insertions(+) create mode 100644 doc/news/changes/minor/20230217tjhei create mode 100644 tests/vector_tools/add_constant_01.cc create mode 100644 tests/vector_tools/add_constant_01.output create mode 100644 tests/vector_tools/add_constant_02.cc create mode 100644 tests/vector_tools/add_constant_02.mpirun=2.output diff --git a/doc/news/changes/minor/20230217tjhei b/doc/news/changes/minor/20230217tjhei new file mode 100644 index 0000000000..9ad1a0bab5 --- /dev/null +++ b/doc/news/changes/minor/20230217tjhei @@ -0,0 +1,4 @@ +New: Implement `VectorTools::add_constant` that adds a constant +to a given component of a finite element solution. +
+(Timo Heister, 2023/02/17) diff --git a/include/deal.II/numerics/vector_tools_mean_value.h b/include/deal.II/numerics/vector_tools_mean_value.h index 2a0646be2c..6e0246832c 100644 --- a/include/deal.II/numerics/vector_tools_mean_value.h +++ b/include/deal.II/numerics/vector_tools_mean_value.h @@ -100,6 +100,38 @@ namespace VectorTools void subtract_mean_value(VectorType &v, const std::vector &p_select = {}); + /** + * Add the constant @p constant_adjustment to the specified + * component of the finite element function given by the coefficient + * vector @p solution defined by the given DoFHandler. + * + * This operation is a common operation to compute a solution with + * mean pressure zero for a Stokes flow problem. Here, one can + * use VectorTools::compute_mean_value() to compute the value to + * subtract. + * + * For a nodal finite element like FE_Q, this function will simply add + * the value @p constant_adjustment to each coefficient of the corresponding + * component. If you have the component in a separate block @p b, you + * could directly use + * solution.block(b) += constant_adjustment + * instead of calling this function (and this would be more efficient). + * For other finite element spaces like FE_DGP, the logic is more + * complicated and handled correctly by this function. + * + * @note Not all kinds of finite elements are supported and the selected + * component must not be part of a non-primitive element for this + * implementation to work. + * + * @note In contrast to subtract_mean_value(), this function can + * adjust a single component of a distributed vector. + */ + template + void + add_constant(VectorType & solution, + const DoFHandler & dof_handler, + const unsigned int component, + const typename VectorType::value_type constant_adjustment); /** * Compute the mean value of one component of the solution. diff --git a/include/deal.II/numerics/vector_tools_mean_value.templates.h b/include/deal.II/numerics/vector_tools_mean_value.templates.h index 5d360f6039..67d451e768 100644 --- a/include/deal.II/numerics/vector_tools_mean_value.templates.h +++ b/include/deal.II/numerics/vector_tools_mean_value.templates.h @@ -18,6 +18,10 @@ #include +#include +#include +#include +#include #include #include @@ -124,6 +128,151 @@ namespace VectorTools } } // namespace internal + + + template + void + add_constant(VectorType & solution, + const DoFHandler & dof_handler, + const unsigned int component, + const typename VectorType::value_type constant_adjustment) + { + Assert(dof_handler.has_hp_capabilities() == false, ExcNotImplemented()); + + AssertDimension(solution.size(), dof_handler.n_dofs()); + AssertIndexRange(component, dof_handler.get_fe().n_components()); + + const FiniteElement &fe_system = dof_handler.get_fe(); + const FiniteElement &fe = fe_system.get_sub_fe(component, 1); + + if ((dynamic_cast *>(&fe) != nullptr)) + { + // The FE to modify is an FE_DGP, which is not a nodal + // element. The first shape function of a DGP element happens + // to be the constant function, so we just have to adjust the + // corresponding DoF on each cell: + std::vector local_dof_indices( + dof_handler.get_fe().dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->get_dof_indices(local_dof_indices); + const unsigned int first_pressure_dof = + fe_system.component_to_system_index(component, 0); + + // Make sure that this DoF is really owned by the + // current processor: + Assert(dof_handler.locally_owned_dofs().is_element( + local_dof_indices[first_pressure_dof]), + ExcInternalError()); + + // Then adjust its value: + solution(local_dof_indices[first_pressure_dof]) += + constant_adjustment; + } + + solution.compress(VectorOperation::add); + } + else if ((dynamic_cast *>(&fe) != nullptr) || + (dynamic_cast *>(&fe) != nullptr)) + { + // We need to make sure to not touch DoFs shared between cells more + // than once. Instead of counting or limiting the number of times + // we touch an individual vector entry, we make a copy of the vector + // and use that as the input and add the constant to it. This way + // it does not matter how often an individual entry is touched. + + VectorType copy(solution); + copy = solution; + + std::vector local_dof_indices( + dof_handler.get_fe().dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->get_dof_indices(local_dof_indices); + for (unsigned i = 0; i < dof_handler.get_fe().dofs_per_cell; ++i) + { + if (!fe_system.is_primitive(i)) + continue; + + const auto component_and_index = + fe_system.system_to_component_index(i); + if (component_and_index.first == component) + { + const types::global_dof_index idx = local_dof_indices[i]; + // Make sure that this DoF is really owned by the + // current processor: + if (!dof_handler.locally_owned_dofs().is_element(idx)) + continue; + + // Then adjust its value: + solution(idx) = copy(idx) + constant_adjustment; + } + } + } + + solution.compress(VectorOperation::insert); + } + else if ((dynamic_cast *>(&fe) != nullptr) || + (dynamic_cast *>(&fe) != + nullptr)) + { + // Add the constant to every single shape function per cell + + std::vector local_dof_indices( + dof_handler.get_fe().dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->get_dof_indices(local_dof_indices); + for (unsigned i = 0; i < dof_handler.get_fe().dofs_per_cell; ++i) + { + if (!fe_system.is_primitive(i)) + continue; + + const auto component_and_index = + fe_system.system_to_component_index(i); + if (component_and_index.first == component) + { + // Make sure that this DoF is really owned by the + // current processor: + Assert(dof_handler.locally_owned_dofs().is_element( + local_dof_indices[i]), + ExcInternalError()); + + // Then adjust its value: + solution(local_dof_indices[i]) += constant_adjustment; + } + } + } + + solution.compress(VectorOperation::add); + } + else + AssertThrow(false, ExcNotImplemented()); + } + + + +#ifdef DEAL_II_WITH_TRILINOS + template + void + add_constant(LinearAlgebra::EpetraWrappers::Vector &, + const DoFHandler &, + const unsigned int, + const double) + { + // TODO: no vector access using operator() + AssertThrow(false, ExcNotImplemented()); + } +#endif + + + template typename VectorType::value_type compute_mean_value( diff --git a/source/numerics/vector_tools_mean_value.inst.in b/source/numerics/vector_tools_mean_value.inst.in index 246a32d601..0085f615a9 100644 --- a/source/numerics/vector_tools_mean_value.inst.in +++ b/source/numerics/vector_tools_mean_value.inst.in @@ -36,6 +36,12 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const VEC &, const unsigned int); + template void + add_constant<>(VEC & solution, + const DoFHandler &dof_handler, + const unsigned int component, + const typename VEC::value_type constant_adjustment); \} #endif } diff --git a/tests/vector_tools/add_constant_01.cc b/tests/vector_tools/add_constant_01.cc new file mode 100644 index 0000000000..8519215e3f --- /dev/null +++ b/tests/vector_tools/add_constant_01.cc @@ -0,0 +1,176 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 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. +// +// --------------------------------------------------------------------- + +// Test VectorTools::add_constant() + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include "../tests.h" + + + +// 0 (dim components) +// 1 +// 2 +// 3 +template +class Reference : public Function +{ +public: + Reference() + : Function(dim + 3) + {} + + double + value(const Point &p, const unsigned int c) const + { + if (c == dim) + return 1.0; + if (c == dim + 1) + return 2.0; + if (c == dim + 2) + return 3.0; + return 0.0; + } +}; + + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FESystem fe(FE_RaviartThomas(1), + 1, + FE_Q(2), + 1, + FE_DGQ(2), + 1, + FE_DGP(1), + 1); + DoFHandler dofhandler(tria); + dofhandler.distribute_dofs(fe); + + VectorType vec(dofhandler.n_dofs()); + vec = 0.0; + for (unsigned int c = dim; c < dim + 3; ++c) + VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1)); + + Vector cellwise_errors(tria.n_active_cells()); + QIterated quadrature(QTrapezoid<1>(), 5); + + const dealii::Function *w = nullptr; + VectorTools::integrate_difference(dofhandler, + vec, + Reference(), + cellwise_errors, + quadrature, + VectorTools::L2_norm); + + const double error = VectorTools::compute_global_error(tria, + cellwise_errors, + VectorTools::L2_norm); + + AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference")); +} + +template +void +test_simplex() +{ + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 4); + + FESystem fe(FE_SimplexP(1), + dim, + FE_SimplexP(2), + 1, + FE_SimplexP(1), + 1, + FE_SimplexDGP(1), + 1); + DoFHandler dofhandler(tria); + dofhandler.distribute_dofs(fe); + + VectorType vec(dofhandler.n_dofs()); + vec = 0.0; + for (unsigned int c = dim; c < dim + 3; ++c) + VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1)); + + Vector cellwise_errors(tria.n_active_cells()); + QGaussSimplex quadrature(2); + + const dealii::Function *w = nullptr; + MappingFE mapping(FE_SimplexP(1)); + VectorTools::integrate_difference(mapping, + dofhandler, + vec, + Reference(), + cellwise_errors, + quadrature, + VectorTools::L2_norm); + + const double error = VectorTools::compute_global_error(tria, + cellwise_errors, + VectorTools::L2_norm); + + AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference")); +} + + +template +void +big() +{ + test>(); + test_simplex>(); +} + + +int +main(int argc, char **argv) +{ + initlog(); + + big<2>(); + big<3>(); + + deallog << "OK" << std::endl; +} diff --git a/tests/vector_tools/add_constant_01.output b/tests/vector_tools/add_constant_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/vector_tools/add_constant_01.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/vector_tools/add_constant_02.cc b/tests/vector_tools/add_constant_02.cc new file mode 100644 index 0000000000..40b2e1178e --- /dev/null +++ b/tests/vector_tools/add_constant_02.cc @@ -0,0 +1,203 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 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. +// +// --------------------------------------------------------------------- + +// Test VectorTools::add_constant() in parallel + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + + + +// 0 (dim components) +// 1 +// 2 +// 3 +template +class Reference : public Function +{ +public: + Reference() + : Function(dim + 3) + {} + + double + value(const Point &p, const unsigned int c) const + { + if (c == dim) + return 1.0; + if (c == dim + 1) + return 2.0; + if (c == dim + 2) + return 3.0; + return 0.0; + } +}; + + + +template +void +test() +{ + parallel::shared::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FESystem fe(FE_RaviartThomas(1), + 1, + FE_Q(2), + 1, + FE_DGQ(2), + 1, + FE_DGP(1), + 1); + DoFHandler dofhandler(tria); + dofhandler.distribute_dofs(fe); + + VectorType vec(dofhandler.locally_owned_dofs(), MPI_COMM_WORLD); + vec = 0.0; + + for (unsigned int c = dim; c < dim + 3; ++c) + VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1)); + + Vector cellwise_errors(tria.n_active_cells()); + QIterated quadrature(QTrapezoid<1>(), 5); + + const dealii::Function *w = nullptr; + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dofhandler, relevant_dofs); + VectorType ghosted(dofhandler.locally_owned_dofs(), + relevant_dofs, + MPI_COMM_WORLD); + ghosted = vec; + VectorTools::integrate_difference(dofhandler, + ghosted, + Reference(), + cellwise_errors, + quadrature, + VectorTools::L2_norm); + + const double error = VectorTools::compute_global_error(tria, + cellwise_errors, + VectorTools::L2_norm); + + AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference")); +} + +template +void +test_simplex() +{ + parallel::shared::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 4); + + FESystem fe(FE_SimplexP(1), + dim, + FE_SimplexP(1), + 1, + FE_SimplexP(1), + 1, + FE_SimplexDGP(1), + 1); + DoFHandler dofhandler(tria); + dofhandler.distribute_dofs(fe); + + VectorType vec(dofhandler.locally_owned_dofs(), MPI_COMM_WORLD); + vec = 0.0; + for (unsigned int c = dim; c < dim + 3; ++c) + VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1)); + + Vector cellwise_errors(tria.n_active_cells()); + QGaussSimplex quadrature(2); + + const dealii::Function *w = nullptr; + MappingFE mapping(FE_SimplexP(1)); + + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dofhandler, relevant_dofs); + VectorType ghosted(dofhandler.locally_owned_dofs(), + relevant_dofs, + MPI_COMM_WORLD); + ghosted = vec; + VectorTools::integrate_difference(mapping, + dofhandler, + ghosted, + Reference(), + cellwise_errors, + quadrature, + VectorTools::L2_norm); + + const double error = VectorTools::compute_global_error(tria, + cellwise_errors, + VectorTools::L2_norm); + + AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference")); +} + + +template +void +big() +{ + test>(); + test_simplex>(); + +#ifdef DEAL_II_WITH_TRILINOS + test(); + test_simplex(); +#endif +#ifdef DEAL_II_WITH_PETSC + test(); + test_simplex(); +#endif +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + + big<2>(); + big<3>(); + + deallog << "OK" << std::endl; +} diff --git a/tests/vector_tools/add_constant_02.mpirun=2.output b/tests/vector_tools/add_constant_02.mpirun=2.output new file mode 100644 index 0000000000..52103938b0 --- /dev/null +++ b/tests/vector_tools/add_constant_02.mpirun=2.output @@ -0,0 +1,5 @@ + +DEAL:0::OK + +DEAL:1::OK + -- 2.39.5