--- /dev/null
+New: Implement `VectorTools::add_constant` that adds a constant
+to a given component of a finite element solution.
+<br>
+(Timo Heister, 2023/02/17)
void
subtract_mean_value(VectorType &v, const std::vector<bool> &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
+ * <code>solution.block(b) += constant_adjustment</code>
+ * 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 <class VectorType, int dim, int spacedim = dim>
+ void
+ add_constant(VectorType & solution,
+ const DoFHandler<dim, spacedim> & dof_handler,
+ const unsigned int component,
+ const typename VectorType::value_type constant_adjustment);
/**
* Compute the mean value of one component of the solution.
#include <deal.II/distributed/tria_base.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/filtered_iterator.h>
}
} // namespace internal
+
+
+ template <class VectorType, int dim, int spacedim>
+ void
+ add_constant(VectorType & solution,
+ const DoFHandler<dim, spacedim> & 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<dim, spacedim> &fe_system = dof_handler.get_fe();
+ const FiniteElement<dim, spacedim> &fe = fe_system.get_sub_fe(component, 1);
+
+ if ((dynamic_cast<const FE_DGP<dim, spacedim> *>(&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<types::global_dof_index> 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<const FE_Q<dim, spacedim> *>(&fe) != nullptr) ||
+ (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&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<types::global_dof_index> 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<const FE_DGQ<dim, spacedim> *>(&fe) != nullptr) ||
+ (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe) !=
+ nullptr))
+ {
+ // Add the constant to every single shape function per cell
+
+ std::vector<types::global_dof_index> 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 <int dim, int spacedim>
+ void
+ add_constant(LinearAlgebra::EpetraWrappers::Vector &,
+ const DoFHandler<dim, spacedim> &,
+ const unsigned int,
+ const double)
+ {
+ // TODO: no vector access using operator()
+ AssertThrow(false, ExcNotImplemented());
+ }
+#endif
+
+
+
template <int dim, typename VectorType, int spacedim>
typename VectorType::value_type
compute_mean_value(
const VEC &,
const unsigned int);
+ template void
+ add_constant<>(VEC & solution,
+ const DoFHandler<deal_II_dimension,
+ deal_II_space_dimension> &dof_handler,
+ const unsigned int component,
+ const typename VEC::value_type constant_adjustment);
\}
#endif
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+// 0 (dim components)
+// 1
+// 2
+// 3
+template <int dim>
+class Reference : public Function<dim>
+{
+public:
+ Reference()
+ : Function<dim>(dim + 3)
+ {}
+
+ double
+ value(const Point<dim> &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 <int dim, class VectorType>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ FESystem<dim> fe(FE_RaviartThomas<dim>(1),
+ 1,
+ FE_Q<dim>(2),
+ 1,
+ FE_DGQ<dim>(2),
+ 1,
+ FE_DGP<dim>(1),
+ 1);
+ DoFHandler<dim> 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<double> cellwise_errors(tria.n_active_cells());
+ QIterated<dim> quadrature(QTrapezoid<1>(), 5);
+
+ const dealii::Function<dim, double> *w = nullptr;
+ VectorTools::integrate_difference(dofhandler,
+ vec,
+ Reference<dim>(),
+ 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 <int dim, class VectorType>
+void
+test_simplex()
+{
+ Triangulation<dim> tria;
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 4);
+
+ FESystem<dim> fe(FE_SimplexP<dim>(1),
+ dim,
+ FE_SimplexP<dim>(2),
+ 1,
+ FE_SimplexP<dim>(1),
+ 1,
+ FE_SimplexDGP<dim>(1),
+ 1);
+ DoFHandler<dim> 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<double> cellwise_errors(tria.n_active_cells());
+ QGaussSimplex<dim> quadrature(2);
+
+ const dealii::Function<dim, double> *w = nullptr;
+ MappingFE<dim> mapping(FE_SimplexP<dim>(1));
+ VectorTools::integrate_difference(mapping,
+ dofhandler,
+ vec,
+ Reference<dim>(),
+ 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 <int dim>
+void
+big()
+{
+ test<dim, Vector<double>>();
+ test_simplex<dim, Vector<double>>();
+}
+
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ big<2>();
+ big<3>();
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+// 0 (dim components)
+// 1
+// 2
+// 3
+template <int dim>
+class Reference : public Function<dim>
+{
+public:
+ Reference()
+ : Function<dim>(dim + 3)
+ {}
+
+ double
+ value(const Point<dim> &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 <int dim, class VectorType>
+void
+test()
+{
+ parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ FESystem<dim> fe(FE_RaviartThomas<dim>(1),
+ 1,
+ FE_Q<dim>(2),
+ 1,
+ FE_DGQ<dim>(2),
+ 1,
+ FE_DGP<dim>(1),
+ 1);
+ DoFHandler<dim> 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<double> cellwise_errors(tria.n_active_cells());
+ QIterated<dim> quadrature(QTrapezoid<1>(), 5);
+
+ const dealii::Function<dim, double> *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<dim>(),
+ 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 <int dim, class VectorType>
+void
+test_simplex()
+{
+ parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 4);
+
+ FESystem<dim> fe(FE_SimplexP<dim>(1),
+ dim,
+ FE_SimplexP<dim>(1),
+ 1,
+ FE_SimplexP<dim>(1),
+ 1,
+ FE_SimplexDGP<dim>(1),
+ 1);
+ DoFHandler<dim> 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<double> cellwise_errors(tria.n_active_cells());
+ QGaussSimplex<dim> quadrature(2);
+
+ const dealii::Function<dim, double> *w = nullptr;
+ MappingFE<dim> mapping(FE_SimplexP<dim>(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<dim>(),
+ 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 <int dim>
+void
+big()
+{
+ test<dim, LinearAlgebra::distributed::Vector<double>>();
+ test_simplex<dim, LinearAlgebra::distributed::Vector<double>>();
+
+#ifdef DEAL_II_WITH_TRILINOS
+ test<dim, TrilinosWrappers::MPI::Vector>();
+ test_simplex<dim, TrilinosWrappers::MPI::Vector>();
+#endif
+#ifdef DEAL_II_WITH_PETSC
+ test<dim, PETScWrappers::MPI::Vector>();
+ test_simplex<dim, PETScWrappers::MPI::Vector>();
+#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;
+}
--- /dev/null
+
+DEAL:0::OK
+
+DEAL:1::OK
+
projects / dealii.git / commitdiff