From e06c1a6748ac71818d4b9900468e19e649e03684 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Fri, 28 Jan 2022 18:25:22 +0100 Subject: [PATCH] Added hp-version for `VectorTools::compute_mean_value()`. --- .../numerics/vector_tools_mean_value.h | 22 +++ .../vector_tools_mean_value.templates.h | 68 ++++++--- ...mean_value.cc => compute_mean_value_01.cc} | 0 ...nos=true.mpirun=10.with_p4est=true.output} | 0 ...inos=true.mpirun=4.with_p4est=true.output} | 0 tests/mpi/compute_mean_value_02.cc | 135 ++++++++++++++++++ ...inos=true.mpirun=10.with_p4est=true.output | 3 + ...linos=true.mpirun=4.with_p4est=true.output | 3 + 8 files changed, 212 insertions(+), 19 deletions(-) rename tests/mpi/{compute_mean_value.cc => compute_mean_value_01.cc} (100%) rename tests/mpi/{compute_mean_value.with_trilinos=true.mpirun=10.with_p4est=true.output => compute_mean_value_01.with_trilinos=true.mpirun=10.with_p4est=true.output} (100%) rename tests/mpi/{compute_mean_value.with_trilinos=true.mpirun=4.with_p4est=true.output => compute_mean_value_01.with_trilinos=true.mpirun=4.with_p4est=true.output} (100%) create mode 100644 tests/mpi/compute_mean_value_02.cc create mode 100644 tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=10.with_p4est=true.output create mode 100644 tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=4.with_p4est=true.output diff --git a/include/deal.II/numerics/vector_tools_mean_value.h b/include/deal.II/numerics/vector_tools_mean_value.h index 15a45479c8..44620e92b9 100644 --- a/include/deal.II/numerics/vector_tools_mean_value.h +++ b/include/deal.II/numerics/vector_tools_mean_value.h @@ -23,8 +23,16 @@ DEAL_II_NAMESPACE_OPEN +#ifndef DOXYGEN +// forward declarations template class DoFHandler; +namespace hp +{ + template + class MappingCollection; +} +#endif namespace VectorTools { @@ -112,6 +120,20 @@ namespace VectorTools */ template typename VectorType::value_type + compute_mean_value( + const hp::MappingCollection &mapping_collection, + const DoFHandler & dof, + const hp::QCollection & q_collection, + const VectorType & v, + const unsigned int component); + + /** + * Calls the other compute_mean_value() function, see above, for the non-hp + * case. That means, it requires a single FiniteElement, a single Quadrature, + * and a single Mapping object. + */ + template + typename VectorType::value_type compute_mean_value(const Mapping & mapping, const DoFHandler &dof, const Quadrature & quadrature, 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 91804a5d35..6266cf4441 100644 --- a/include/deal.II/numerics/vector_tools_mean_value.templates.h +++ b/include/deal.II/numerics/vector_tools_mean_value.templates.h @@ -21,6 +21,10 @@ #include +#include +#include +#include + #include #include #include @@ -123,24 +127,29 @@ namespace VectorTools template typename VectorType::value_type - compute_mean_value(const Mapping & mapping, - const DoFHandler &dof, - const Quadrature & quadrature, - const VectorType & v, - const unsigned int component) + compute_mean_value( + const hp::MappingCollection &mapping_collection, + const DoFHandler & dof, + const hp::QCollection & q_collection, + const VectorType & v, + const unsigned int component) { using Number = typename VectorType::value_type; - Assert(v.size() == dof.n_dofs(), - ExcDimensionMismatch(v.size(), dof.n_dofs())); - AssertIndexRange(component, dof.get_fe(0).n_components()); - FEValues fe(mapping, - dof.get_fe(), - quadrature, - UpdateFlags(update_JxW_values | update_values)); + const hp::FECollection &fe_collection = + dof.get_fe_collection(); + const unsigned int n_components = fe_collection.n_components(); - std::vector> values( - quadrature.size(), Vector(dof.get_fe(0).n_components())); + AssertDimension(v.size(), dof.n_dofs()); + AssertIndexRange(component, n_components); + + hp::FEValues fe_values_collection( + mapping_collection, + fe_collection, + q_collection, + UpdateFlags(update_JxW_values | update_values)); + + std::vector> values; Number mean = Number(); typename numbers::NumberTraits::real_type area = 0.; @@ -148,12 +157,17 @@ namespace VectorTools for (const auto &cell : dof.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) { - fe.reinit(cell); - fe.get_function_values(v, values); - for (unsigned int k = 0; k < quadrature.size(); ++k) + fe_values_collection.reinit(cell); + const FEValues &fe_values = + fe_values_collection.get_present_fe_values(); + + values.resize(fe_values.n_quadrature_points, + Vector(n_components)); + fe_values.get_function_values(v, values); + for (unsigned int k = 0; k < fe_values.n_quadrature_points; ++k) { - mean += fe.JxW(k) * values[k](component); - area += fe.JxW(k); + mean += fe_values.JxW(k) * values[k](component); + area += fe_values.JxW(k); } } @@ -192,6 +206,22 @@ namespace VectorTools } + template + typename VectorType::value_type + compute_mean_value(const Mapping & mapping, + const DoFHandler &dof, + const Quadrature & quadrature, + const VectorType & v, + const unsigned int component) + { + return compute_mean_value(hp::MappingCollection(mapping), + dof, + hp::QCollection(quadrature), + v, + component); + } + + template typename VectorType::value_type compute_mean_value(const DoFHandler &dof, diff --git a/tests/mpi/compute_mean_value.cc b/tests/mpi/compute_mean_value_01.cc similarity index 100% rename from tests/mpi/compute_mean_value.cc rename to tests/mpi/compute_mean_value_01.cc diff --git a/tests/mpi/compute_mean_value.with_trilinos=true.mpirun=10.with_p4est=true.output b/tests/mpi/compute_mean_value_01.with_trilinos=true.mpirun=10.with_p4est=true.output similarity index 100% rename from tests/mpi/compute_mean_value.with_trilinos=true.mpirun=10.with_p4est=true.output rename to tests/mpi/compute_mean_value_01.with_trilinos=true.mpirun=10.with_p4est=true.output diff --git a/tests/mpi/compute_mean_value.with_trilinos=true.mpirun=4.with_p4est=true.output b/tests/mpi/compute_mean_value_01.with_trilinos=true.mpirun=4.with_p4est=true.output similarity index 100% rename from tests/mpi/compute_mean_value.with_trilinos=true.mpirun=4.with_p4est=true.output rename to tests/mpi/compute_mean_value_01.with_trilinos=true.mpirun=4.with_p4est=true.output diff --git a/tests/mpi/compute_mean_value_02.cc b/tests/mpi/compute_mean_value_02.cc new file mode 100644 index 0000000000..52bda4e4f3 --- /dev/null +++ b/tests/mpi/compute_mean_value_02.cc @@ -0,0 +1,135 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 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 VectorTools::compute_mean_value for parallel computations for +// hp-adaptive applications. +// +// we interpolate a linear function onto the grid with a symmetric mesh. +// the mean value of the interpolation must be the mean of the linear function + +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + + + +template +class LinearFunction : public Function +{ +public: + double + value(const Point &p, const unsigned int) const + { + return p[0] + 2; + } +}; + + +template +void +test() +{ + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + { + GridGenerator::hyper_ball(tr); + tr.refine_global(2); + } + + hp::MappingCollection mapping_collection((MappingQ1())); + + hp::FECollection fe_collection; + hp::QCollection q_collection; + for (unsigned int degree = 1; degree <= 3; ++degree) + { + fe_collection.push_back(FE_Q(degree)); + q_collection.push_back(QGauss(degree + 1)); + } + + DoFHandler dofh(tr); + { + // set fe indices arbitrarily + unsigned int i = 0; + for (const auto &cell : + dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + cell->set_active_fe_index(i++ % fe_collection.size()); + dofh.distribute_dofs(fe_collection); + } + + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs(dofh, relevant_set); + TrilinosWrappers::MPI::Vector x_rel(relevant_set, dofh.get_communicator()); + { + TrilinosWrappers::MPI::Vector interpolated(dofh.locally_owned_dofs(), + dofh.get_communicator()); + VectorTools::interpolate(dofh, LinearFunction(), interpolated); + x_rel = interpolated; + } + + const double mean = VectorTools::compute_mean_value( + mapping_collection, dofh, q_collection, x_rel, 0); + + deallog << "mean=" << mean << std::endl; + + Assert(std::fabs(mean - 2) < 1e-3, ExcInternalError()); +} + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=10.with_p4est=true.output b/tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=10.with_p4est=true.output new file mode 100644 index 0000000000..75c95df4c9 --- /dev/null +++ b/tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=10.with_p4est=true.output @@ -0,0 +1,3 @@ + +DEAL:0:2d::mean=2.00000 +DEAL:0:3d::mean=2.00000 diff --git a/tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=4.with_p4est=true.output b/tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=4.with_p4est=true.output new file mode 100644 index 0000000000..75c95df4c9 --- /dev/null +++ b/tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=4.with_p4est=true.output @@ -0,0 +1,3 @@ + +DEAL:0:2d::mean=2.00000 +DEAL:0:3d::mean=2.00000 -- 2.39.5