]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow using VectorTools::subtract_mean_value for all distributed vector classes 4463/head
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Thu, 1 Jun 2017 15:57:22 +0000 (17:57 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Thu, 1 Jun 2017 21:41:36 +0000 (23:41 +0200)
include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/vector_tools_mean_value.cc
tests/numerics/subtract_mean_value_01.cc [new file with mode: 0644]
tests/numerics/subtract_mean_value_01.output [new file with mode: 0644]
tests/numerics/subtract_mean_value_02.cc [new file with mode: 0644]
tests/numerics/subtract_mean_value_02.with_mpi=true.mpirun=2.output [new file with mode: 0644]
tests/petsc/subtract_mean_value_03.cc [new file with mode: 0644]
tests/petsc/subtract_mean_value_03.with_petsc=true.with_mpi=true.mpirun=2.output [new file with mode: 0644]
tests/trilinos/subtract_mean_value_04.cc [new file with mode: 0644]
tests/trilinos/subtract_mean_value_04.with_trilinos=true.with_mpi=true.mpirun=2.output [new file with mode: 0644]

index b4c06996025af5d89d6f649105fc4f3ab89dd43f..bdd53a4df144d20f8fdd88069fc718e773ecb38d 100644 (file)
@@ -2830,8 +2830,8 @@ namespace VectorTools
    * not equal to $(1,1,\ldots,1)^T$. For such elements, a different procedure
    * has to be used when subtracting the mean value.
    *
-   * @warning This function is only implemented for Vector and BlockVector. It
-   * is not implemented for any of the distributed vector classes.
+   * @warning This function can only be used for distributed vector classes
+   * provided the boolean mask is empty, i.e. selecting the whole vector.
    */
   template <typename VectorType>
   void subtract_mean_value(VectorType              &v,
index 7ee26717ae5eeac9fd888d9f0b4e5236feae0784..482940d7697aa5a53087085ee33179ecece14e82 100644 (file)
@@ -7753,48 +7753,67 @@ namespace VectorTools
     return gradient[0];
   }
 
+  namespace
+  {
+    template <typename VectorType>
+    typename std::enable_if<dealii::is_serial_vector< VectorType >::value == true>::type
+    internal_subtract_mean_value(VectorType              &v,
+                                 const std::vector<bool> &p_select)
+    {
+      if (p_select.size() == 0)
+        {
+          // In case of an empty boolean mask operate on the whole vector:
+          v.add( - v.mean_value() );
+        }
+      else
+        {
+          const unsigned int n = v.size();
 
+          Assert(p_select.size() == n,
+                 ExcDimensionMismatch(p_select.size(), n));
 
-  template <typename VectorType>
-  void
-  subtract_mean_value(VectorType              &v,
-                      const std::vector<bool> &p_select)
-  {
-    if (p_select.size() == 0)
-      {
-        // In case of an empty boolean mask operate on the whole vector:
-        v.add( - v.mean_value() );
-      }
-    else
-      {
-        // This function is not implemented for distributed vectors.
-        Assert((dealii::is_serial_vector< VectorType >::value == true),
-               ExcNotImplemented());
+          typename VectorType::value_type s = 0.;
+          unsigned int counter = 0;
+          for (unsigned int i=0; i<n; ++i)
+            if (p_select[i])
+              {
+                typename VectorType::value_type vi = v(i);
+                s += vi;
+                ++counter;
+              }
+          // Error out if we have not constrained anything. Note that in this
+          // case the vector v is always nonempty.
+          Assert (n == 0 || counter > 0, ComponentMask::ExcNoComponentSelected());
 
-        const unsigned int n = v.size();
+          s /= counter;
 
-        Assert(p_select.size() == n,
-               ExcDimensionMismatch(p_select.size(), n));
+          for (unsigned int i=0; i<n; ++i)
+            if (p_select[i])
+              v(i) -= s;
+        }
+    }
 
-        typename VectorType::value_type s = 0.;
-        unsigned int counter = 0;
-        for (unsigned int i=0; i<n; ++i)
-          if (p_select[i])
-            {
-              typename VectorType::value_type vi = v(i);
-              s += vi;
-              ++counter;
-            }
-        // Error out if we have not constrained anything. Note that in this
-        // case the vector v is always nonempty.
-        Assert (n == 0 || counter > 0, ComponentMask::ExcNoComponentSelected());
 
-        s /= counter;
 
-        for (unsigned int i=0; i<n; ++i)
-          if (p_select[i])
-            v(i) -= s;
-      }
+    template <typename VectorType>
+    typename std::enable_if<dealii::is_serial_vector< VectorType >::value == false>::type
+    internal_subtract_mean_value(VectorType              &v,
+                                 const std::vector<bool> &p_select)
+    {
+      (void) p_select;
+      Assert(p_select.size() == 0, ExcNotImplemented());
+      // In case of an empty boolean mask operate on the whole vector:
+      v.add( - v.mean_value() );
+    }
+  }
+
+
+  template <typename VectorType>
+  void
+  subtract_mean_value(VectorType              &v,
+                      const std::vector<bool> &p_select)
+  {
+    internal_subtract_mean_value(v,p_select);
   }
 
   namespace
index 2861b76a98ef8f1867a6b47148c18869a39c3870..43b8e51e983a977818377df0668a21e091c6b670 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
-#if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
-namespace VectorTools
-{
-  template <>
-  void
-  subtract_mean_value(LinearAlgebra::EpetraWrappers::Vector &,
-                      const std::vector<bool> &)
-  {
-    Assert(false,ExcNotImplemented());
-  }
-}
-#endif
-
 // ---------------------------- explicit instantiations --------------------
 #include "vector_tools_mean_value.inst"
 
diff --git a/tests/numerics/subtract_mean_value_01.cc b/tests/numerics/subtract_mean_value_01.cc
new file mode 100644 (file)
index 0000000..3a938d9
--- /dev/null
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check VectorTools::subtract_mean_value() for deal.II serial vectors
+
+#include "../tests.h"
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/la_vector.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+template <class VectorType>
+void test (VectorType &v)
+{
+  std::vector<bool> filter(v.size(), false);
+  // set some elements of the vector
+  for (unsigned int i=0; i<v.size(); i+=1+i)
+    {
+      filter[i] = true;
+      v(i) = i;
+    }
+
+  // then check the norm
+  VectorTools::subtract_mean_value(v, filter);
+  AssertThrow (std::fabs(v.mean_value()) < 1e-10*v.l2_norm(),
+               ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  initlog();
+
+  try
+    {
+      {
+        Vector<double> v (10);
+        test (v);
+      }
+
+      {
+        Vector<float>  v (10);
+        test (v);
+      }
+
+      {
+        BlockVector<double> v (std::vector<types::global_dof_index>(1,10));
+        test (v);
+      }
+
+      {
+        BlockVector<float> v (std::vector< types::global_dof_index>(1,10));
+        test (v);
+      }
+
+      {
+        LinearAlgebra::Vector<double> v(10);
+        test (v);
+      }
+
+      {
+        LinearAlgebra::Vector<float>  v(10);
+        test (v);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/numerics/subtract_mean_value_01.output b/tests/numerics/subtract_mean_value_01.output
new file mode 100644 (file)
index 0000000..0aa61ff
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/numerics/subtract_mean_value_02.cc b/tests/numerics/subtract_mean_value_02.cc
new file mode 100644 (file)
index 0000000..828181b
--- /dev/null
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check VectorTools::subtract_mean_value() for deal.II parallel vector
+
+#include "../tests.h"
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+template <class VectorType>
+void test (VectorType &v)
+{
+  // set some elements of the vector
+  unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  for (unsigned int i=5*my_id; i<5*(my_id+1); ++i)
+    {
+      v(i) = i;
+    }
+  v.compress (VectorOperation::insert);
+
+  // then check the norm
+  VectorTools::subtract_mean_value(v);
+  AssertThrow (std::fabs(v.mean_value()) < 1e-10*v.l2_norm(),
+               ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  mpi_initlog();
+
+  unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  try
+    {
+      IndexSet local_range(10);
+      local_range.add_range(5*my_id,5*(my_id+1));
+      IndexSet ghost_indices(10);
+      ghost_indices.add_range(3, 8);
+      {
+        LinearAlgebra::distributed::Vector<double> v
+        (local_range, ghost_indices, MPI_COMM_WORLD);
+        test (v);
+      }
+
+      {
+        LinearAlgebra::distributed::Vector<float>  v
+        (local_range, ghost_indices, MPI_COMM_WORLD);
+        test (v);
+      }
+
+      {
+        LinearAlgebra::distributed::BlockVector<double> v
+        (std::vector<IndexSet>(1, local_range),
+         std::vector<IndexSet>(1, ghost_indices), MPI_COMM_WORLD);
+        test (v);
+      }
+
+      {
+        LinearAlgebra::distributed::BlockVector<float> v
+        (std::vector<IndexSet>(1, local_range),
+         std::vector<IndexSet>(1, ghost_indices), MPI_COMM_WORLD);
+        test (v);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/numerics/subtract_mean_value_02.with_mpi=true.mpirun=2.output b/tests/numerics/subtract_mean_value_02.with_mpi=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..5f42bb2
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/petsc/subtract_mean_value_03.cc b/tests/petsc/subtract_mean_value_03.cc
new file mode 100644 (file)
index 0000000..5c76824
--- /dev/null
@@ -0,0 +1,94 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check VectorTools::subtract_mean_value() for PETSc vectors
+
+#include "../tests.h"
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+template <class VectorType>
+void test (VectorType &v)
+{
+  // set some elements of the vector
+  unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  for (unsigned int i=5*my_id; i<5*(my_id+1); ++i)
+    {
+      v(i) = i;
+    }
+  v.compress (VectorOperation::insert);
+
+  // then check the norm
+  VectorTools::subtract_mean_value(v);
+  AssertThrow (std::fabs(v.mean_value()) < 1e-10*v.l2_norm(),
+               ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+int main (int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  mpi_initlog();
+
+  unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  try
+    {
+      IndexSet local_range(10);
+      local_range.add_range(5*my_id,5*(my_id+1));
+      {
+        PETScWrappers::MPI::Vector v
+        (local_range, MPI_COMM_WORLD);
+        test (v);
+      }
+
+      {
+        PETScWrappers::MPI::BlockVector v
+        (std::vector<IndexSet>(1, local_range), MPI_COMM_WORLD);
+        test (v);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/petsc/subtract_mean_value_03.with_petsc=true.with_mpi=true.mpirun=2.output b/tests/petsc/subtract_mean_value_03.with_petsc=true.with_mpi=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..8b3b075
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::OK
diff --git a/tests/trilinos/subtract_mean_value_04.cc b/tests/trilinos/subtract_mean_value_04.cc
new file mode 100644 (file)
index 0000000..f52250b
--- /dev/null
@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check VectorTools::subtract_mean_value() for Trilinos vectors
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+template <class VectorType>
+void test (VectorType &v)
+{
+  // set some elements of the vector
+  unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  for (unsigned int i=5*my_id; i<5*(my_id+1); ++i)
+    {
+      v(i) = i;
+    }
+  v.compress (VectorOperation::insert);
+
+  // then check the norm
+  VectorTools::subtract_mean_value(v);
+  AssertThrow (std::fabs(v.mean_value()) < 1e-10*v.l2_norm(),
+               ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+int main (int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  mpi_initlog();
+
+  unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  try
+    {
+      IndexSet local_range(10);
+      local_range.add_range(5*my_id,5*(my_id+1));
+      {
+        TrilinosWrappers::MPI::Vector v1 (local_range, MPI_COMM_WORLD);
+        test(v1);
+
+        LinearAlgebra::ReadWriteVector<double> v_tmp (local_range);
+        LinearAlgebra::EpetraWrappers::Vector v2 (local_range, MPI_COMM_WORLD);
+        v_tmp.import(v1, VectorOperation::insert);
+        v2.import(v_tmp, VectorOperation::insert);
+        VectorTools::subtract_mean_value(v2);
+        AssertThrow (std::fabs(v2.mean_value()) < 1e-10*v2.l2_norm(),
+                     ExcInternalError());
+        deallog << "OK" << std::endl;
+      }
+      {
+        TrilinosWrappers::MPI::BlockVector v (std::vector<IndexSet>(1, local_range),
+                                              MPI_COMM_WORLD);
+        test(v);
+      }
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/subtract_mean_value_04.with_trilinos=true.with_mpi=true.mpirun=2.output b/tests/trilinos/subtract_mean_value_04.with_trilinos=true.with_mpi=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..fb71de2
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK

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.