* 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,
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
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"
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
projects / dealii.git / commitdiff