]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added hp-version for `VectorTools::compute_mean_value()`. 13309/head
authorMarc Fehling <mafehling.git@gmail.com>
Fri, 28 Jan 2022 17:25:22 +0000 (18:25 +0100)
committerMarc Fehling <mafehling.git@gmail.com>
Sat, 29 Jan 2022 16:19:51 +0000 (09:19 -0700)
include/deal.II/numerics/vector_tools_mean_value.h
include/deal.II/numerics/vector_tools_mean_value.templates.h
tests/mpi/compute_mean_value_01.cc [moved from tests/mpi/compute_mean_value.cc with 100% similarity]
tests/mpi/compute_mean_value_01.with_trilinos=true.mpirun=10.with_p4est=true.output [moved from tests/mpi/compute_mean_value.with_trilinos=true.mpirun=10.with_p4est=true.output with 100% similarity]
tests/mpi/compute_mean_value_01.with_trilinos=true.mpirun=4.with_p4est=true.output [moved from tests/mpi/compute_mean_value.with_trilinos=true.mpirun=4.with_p4est=true.output with 100% similarity]
tests/mpi/compute_mean_value_02.cc [new file with mode: 0644]
tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=10.with_p4est=true.output [new file with mode: 0644]
tests/mpi/compute_mean_value_02.with_trilinos=true.mpirun=4.with_p4est=true.output [new file with mode: 0644]

index 15a45479c81daa8609f62ab3759709c75aff63ea..44620e92b9c6011039223e243ca870896aecc20c 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+#ifndef DOXYGEN
+// forward declarations
 template <int dim, int spacedim>
 class DoFHandler;
+namespace hp
+{
+  template <int dim, int spacedim>
+  class MappingCollection;
+}
+#endif
 
 namespace VectorTools
 {
@@ -112,6 +120,20 @@ namespace VectorTools
    */
   template <int dim, typename VectorType, int spacedim>
   typename VectorType::value_type
+  compute_mean_value(
+    const hp::MappingCollection<dim, spacedim> &mapping_collection,
+    const DoFHandler<dim, spacedim> &           dof,
+    const hp::QCollection<dim> &                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 <int dim, typename VectorType, int spacedim>
+  typename VectorType::value_type
   compute_mean_value(const Mapping<dim, spacedim> &   mapping,
                      const DoFHandler<dim, spacedim> &dof,
                      const Quadrature<dim> &          quadrature,
index 91804a5d35d8a590f78f6563c0d5c2e3ca969e91..6266cf44418fc2e5854cf5310f54384f674b39dd 100644 (file)
 
 #include <deal.II/grid/filtered_iterator.h>
 
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/q_collection.h>
+
 #include <deal.II/lac/block_vector.h>
 #include <deal.II/lac/la_parallel_block_vector.h>
 #include <deal.II/lac/la_parallel_vector.h>
@@ -123,24 +127,29 @@ namespace VectorTools
 
   template <int dim, typename VectorType, int spacedim>
   typename VectorType::value_type
-  compute_mean_value(const Mapping<dim, spacedim> &   mapping,
-                     const DoFHandler<dim, spacedim> &dof,
-                     const Quadrature<dim> &          quadrature,
-                     const VectorType &               v,
-                     const unsigned int               component)
+  compute_mean_value(
+    const hp::MappingCollection<dim, spacedim> &mapping_collection,
+    const DoFHandler<dim, spacedim> &           dof,
+    const hp::QCollection<dim> &                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<dim, spacedim> fe(mapping,
-                               dof.get_fe(),
-                               quadrature,
-                               UpdateFlags(update_JxW_values | update_values));
+    const hp::FECollection<dim, spacedim> &fe_collection =
+      dof.get_fe_collection();
+    const unsigned int n_components = fe_collection.n_components();
 
-    std::vector<Vector<Number>> values(
-      quadrature.size(), Vector<Number>(dof.get_fe(0).n_components()));
+    AssertDimension(v.size(), dof.n_dofs());
+    AssertIndexRange(component, n_components);
+
+    hp::FEValues<dim, spacedim> fe_values_collection(
+      mapping_collection,
+      fe_collection,
+      q_collection,
+      UpdateFlags(update_JxW_values | update_values));
+
+    std::vector<Vector<Number>> values;
 
     Number                                            mean = Number();
     typename numbers::NumberTraits<Number>::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<dim, spacedim> &fe_values =
+          fe_values_collection.get_present_fe_values();
+
+        values.resize(fe_values.n_quadrature_points,
+                      Vector<Number>(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 <int dim, typename VectorType, int spacedim>
+  typename VectorType::value_type
+  compute_mean_value(const Mapping<dim, spacedim> &   mapping,
+                     const DoFHandler<dim, spacedim> &dof,
+                     const Quadrature<dim> &          quadrature,
+                     const VectorType &               v,
+                     const unsigned int               component)
+  {
+    return compute_mean_value(hp::MappingCollection<dim, spacedim>(mapping),
+                              dof,
+                              hp::QCollection<dim>(quadrature),
+                              v,
+                              component);
+  }
+
+
   template <int dim, typename VectorType, int spacedim>
   typename VectorType::value_type
   compute_mean_value(const DoFHandler<dim, spacedim> &dof,
diff --git a/tests/mpi/compute_mean_value_02.cc b/tests/mpi/compute_mean_value_02.cc
new file mode 100644 (file)
index 0000000..52bda4e
--- /dev/null
@@ -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 <deal.II/base/function.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class LinearFunction : public Function<dim>
+{
+public:
+  double
+  value(const Point<dim> &p, const unsigned int) const
+  {
+    return p[0] + 2;
+  }
+};
+
+
+template <int dim>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
+  {
+    GridGenerator::hyper_ball(tr);
+    tr.refine_global(2);
+  }
+
+  hp::MappingCollection<dim> mapping_collection((MappingQ1<dim>()));
+
+  hp::FECollection<dim> fe_collection;
+  hp::QCollection<dim>  q_collection;
+  for (unsigned int degree = 1; degree <= 3; ++degree)
+    {
+      fe_collection.push_back(FE_Q<dim>(degree));
+      q_collection.push_back(QGauss<dim>(degree + 1));
+    }
+
+  DoFHandler<dim> 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<dim>(), 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 (file)
index 0000000..75c95df
--- /dev/null
@@ -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 (file)
index 0000000..75c95df
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:0:2d::mean=2.00000
+DEAL:0:3d::mean=2.00000

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.