const DataAccessOnFaces src_vector_face_access =
DataAccessOnFaces::unspecified) const;
+ /**
+ * This function is similar to the loop method above, but adds two additional
+ * functors to execute some additional work before and after the cell, face
+ * and boundary integrals are computed.
+ *
+ * The two additional functors work on a range of degrees of freedom,
+ * expressed in terms of the degree-of-freedom numbering of the selected
+ * DoFHandler `dof_handler_index_pre_post` in MPI-local indices. The
+ * arguments to the functors represent a range of degrees of freedom at a
+ * granularity of
+ * internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector entries
+ * (except for the last chunk which is set to the number of locally owned
+ * entries) in the form `[first, last)`. The idea of these functors is to
+ * bring operations on vectors closer to the point where they accessed in a
+ * matrix-free loop, with the goal to increase cache hits by temporal
+ * locality. This loop guarantees that the `operation_before_loop` hits all
+ * relevant unknowns before they are first touched by any of the cell, face or
+ * boundary operations (including the MPI data exchange), allowing to execute
+ * some vector update that the `src` vector depends upon. The
+ * `operation_after_loop` is similar - it starts to execute on a range of DoFs
+ * once all DoFs in that range have been touched for the last time by the
+ * cell, face and boundary operations (including the MPI data exchange),
+ * allowing e.g. to compute some vector operations that depend on the result
+ * of the current cell loop in `dst` or want to modify `src`. The efficiency
+ * of caching depends on the numbering of the degrees of freedom because of
+ * the granularity of the ranges.
+ *
+ * @param cell_operation Pointer to member function of `CLASS` with the
+ * signature <tt>cell_operation (const MatrixFree<dim,Number> &, OutVector &,
+ * InVector &, std::pair<unsigned int,unsigned int> &)</tt> where the first
+ * argument passes the data of the calling class and the last argument
+ * defines the range of cells which should be worked on (typically more than
+ * one cell should be worked on in order to reduce overheads).
+ *
+ * @param face_operation Pointer to member function of `CLASS` with the
+ * signature <tt>face_operation (const MatrixFree<dim,Number> &, OutVector &,
+ * InVector &, std::pair<unsigned int,unsigned int> &)</tt> in analogy to
+ * `cell_operation`, but now the part associated to the work on interior
+ * faces. Note that the MatrixFree framework treats periodic faces as
+ * interior ones, so they will be assigned their correct neighbor after
+ * applying periodicity constraints within the face_operation calls.
+ *
+ * @param boundary_operation Pointer to member function of `CLASS` with the
+ * signature <tt>boundary_operation (const MatrixFree<dim,Number> &, OutVector
+ * &, InVector &, std::pair<unsigned int,unsigned int> &)</tt> in analogy to
+ * `cell_operation` and `face_operation`, but now the part associated to the
+ * work on boundary faces. Boundary faces are separated by their
+ * `boundary_id` and it is possible to query that id using
+ * MatrixFree::get_boundary_id(). Note that both interior and faces use the
+ * same numbering, and faces in the interior are assigned lower numbers than
+ * the boundary faces.
+ *
+ * @param owning_class The object which provides the `cell_operation`
+ * call. To be compatible with this interface, the class must allow to call
+ * `owning_class->cell_operation(...)`.
+ *
+ * @param dst Destination vector holding the result. If the vector is of
+ * type LinearAlgebra::distributed::Vector (or composite objects thereof
+ * such as LinearAlgebra::distributed::BlockVector), the loop calls
+ * LinearAlgebra::distributed::Vector::compress() at the end of the call
+ * internally. For other vectors, including parallel Trilinos or PETSc
+ * vectors, no such call is issued. Note that Trilinos/Epetra or PETSc
+ * vectors do currently not work in parallel because the present class uses
+ * MPI-local index addressing, as opposed to the global addressing implied
+ * by those external libraries.
+ *
+ * @param src Input vector. If the vector is of type
+ * LinearAlgebra::distributed::Vector (or composite objects thereof such as
+ * LinearAlgebra::distributed::BlockVector), the loop calls
+ * LinearAlgebra::distributed::Vector::update_ghost_values() at the start of
+ * the call internally to make sure all necessary data is locally
+ * available. Note, however, that the vector is reset to its original state
+ * at the end of the loop, i.e., if the vector was not ghosted upon entry of
+ * the loop, it will not be ghosted upon finishing the loop.
+ *
+ * @param operation_before_loop This functor can be used to perform an
+ * operation on entries of the `src` and `dst` vectors (or other vectors)
+ * before the operation on cells first touches a particular DoF according to
+ * the general description in the text above. This function is passed a
+ * range of the locally owned degrees of freedom on the selected
+ * `dof_handler_index_pre_post` (in MPI-local numbering).
+ *
+ * @param operation_after_loop This functor can be used to perform an
+ * operation on entries of the `src` and `dst` vectors (or other vectors)
+ * after the operation on cells last touches a particular DoF according to
+ * the general description in the text above. This function is passed a
+ * range of the locally owned degrees of freedom on the selected
+ * `dof_handler_index_pre_post` (in MPI-local numbering).
+ *
+ * @param dof_handler_index_pre_post Since MatrixFree can be initialized
+ * with a vector of DoFHandler objects, each of them will in general have
+ * different vector sizes and thus different ranges returned to
+ * `operation_before_loop` and `operation_after_loop`. Use this variable to
+ * specify which one of the DoFHandler objects the index range should be
+ * associated to. Defaults to the `dof_handler_index` 0.
+ *
+ * @param dst_vector_face_access Set the type of access into the vector
+ * `dst` that will happen inside the body of the @p face_operation
+ * function. As explained in the description of the DataAccessOnFaces
+ * struct, the purpose of this selection is to reduce the amount of data
+ * that must be exchanged over the MPI network (or via `memcpy` if within
+ * the shared memory region of a node) to gain performance. Note that there
+ * is no way to communicate this setting with the FEFaceEvaluation class,
+ * therefore this selection must be made at this site in addition to what is
+ * implemented inside the `face_operation` function. As a consequence, there
+ * is also no way to check that the setting passed to this call is
+ * consistent with what is later done by `FEFaceEvaluation`, and it is the
+ * user's responsibility to ensure correctness of data.
+ *
+ * @param src_vector_face_access Set the type of access into the vector
+ * `src` that will happen inside the body of the @p face_operation function,
+ * in analogy to `dst_vector_face_access`.
+ *
+ * @note The close locality of the `operation_before_loop` and
+ * `operation_after_loop` is currently only implemented for the MPI-only
+ * case. In case threading is enabled, the complete `operation_before_loop`
+ * is scheduled before the parallel loop, and `operation_after_loop` is
+ * scheduled strictly afterwards, due to the complicated dependencies.
+ */
+ template <typename CLASS, typename OutVector, typename InVector>
+ void
+ loop(
+ void (CLASS::*cell_operation)(const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)
+ const,
+ void (CLASS::*face_operation)(const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)
+ const,
+ void (CLASS::*boundary_operation)(
+ const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &) const,
+ const CLASS * owning_class,
+ OutVector & dst,
+ const InVector &src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_loop,
+ const std::function<void(const unsigned int, const unsigned int)>
+ & operation_after_loop,
+ const unsigned int dof_handler_index_pre_post = 0,
+ const DataAccessOnFaces dst_vector_face_access =
+ DataAccessOnFaces::unspecified,
+ const DataAccessOnFaces src_vector_face_access =
+ DataAccessOnFaces::unspecified) const;
+
+ /**
+ * Same as above, but for class member functions which are non-const.
+ */
+ template <typename CLASS, typename OutVector, typename InVector>
+ void
+ loop(void (CLASS::*cell_operation)(
+ const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &),
+ void (CLASS::*face_operation)(
+ const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &),
+ void (CLASS::*boundary_operation)(
+ const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &),
+ const CLASS * owning_class,
+ OutVector & dst,
+ const InVector &src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_loop,
+ const std::function<void(const unsigned int, const unsigned int)>
+ & operation_after_loop,
+ const unsigned int dof_handler_index_pre_post = 0,
+ const DataAccessOnFaces dst_vector_face_access =
+ DataAccessOnFaces::unspecified,
+ const DataAccessOnFaces src_vector_face_access =
+ DataAccessOnFaces::unspecified) const;
+
+ /**
+ * Same as above, but taking an `std::function` as the `cell_operation`,
+ * `face_operation` and `boundary_operation` rather than a class member
+ * function.
+ */
+ template <typename OutVector, typename InVector>
+ void
+ loop(const std::function<
+ void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)> &cell_operation,
+ const std::function<
+ void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)> &face_operation,
+ const std::function<void(
+ const MatrixFree<dim, Number, VectorizedArrayType> &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
+ OutVector & dst,
+ const InVector & src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_loop,
+ const std::function<void(const unsigned int, const unsigned int)>
+ & operation_after_loop,
+ const unsigned int dof_handler_index_pre_post = 0,
+ const DataAccessOnFaces dst_vector_face_access =
+ DataAccessOnFaces::unspecified,
+ const DataAccessOnFaces src_vector_face_access =
+ DataAccessOnFaces::unspecified) const;
+
/**
* This method runs the loop over all cells (in parallel) similarly as
* cell_loop() does. However, this function is intended to be used
+template <int dim, typename Number, typename VectorizedArrayType>
+template <typename OutVector, typename InVector>
+inline void
+MatrixFree<dim, Number, VectorizedArrayType>::loop(
+ const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)>
+ &cell_operation,
+ const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)>
+ &face_operation,
+ const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)>
+ & boundary_operation,
+ OutVector & dst,
+ const InVector &src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_loop,
+ const std::function<void(const unsigned int, const unsigned int)>
+ & operation_after_loop,
+ const unsigned int dof_handler_index_pre_post,
+ const DataAccessOnFaces dst_vector_face_access,
+ const DataAccessOnFaces src_vector_face_access) const
+{
+ using Wrapper =
+ internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
+ InVector,
+ OutVector>;
+ Wrapper wrap(cell_operation, face_operation, boundary_operation);
+ internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
+ InVector,
+ OutVector,
+ Wrapper,
+ true>
+ worker(*this,
+ src,
+ dst,
+ false,
+ wrap,
+ &Wrapper::cell_integrator,
+ &Wrapper::face_integrator,
+ &Wrapper::boundary_integrator,
+ src_vector_face_access,
+ dst_vector_face_access,
+ operation_before_loop,
+ operation_after_loop,
+ dof_handler_index_pre_post);
+
+ task_info.loop(worker);
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+template <typename CLASS, typename OutVector, typename InVector>
+inline void
+MatrixFree<dim, Number, VectorizedArrayType>::loop(
+ void (CLASS::*cell_operation)(const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)
+ const,
+ void (CLASS::*face_operation)(const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &)
+ const,
+ void (CLASS::*boundary_operation)(
+ const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &) const,
+ const CLASS * owning_class,
+ OutVector & dst,
+ const InVector &src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_loop,
+ const std::function<void(const unsigned int, const unsigned int)>
+ & operation_after_loop,
+ const unsigned int dof_handler_index_pre_post,
+ const DataAccessOnFaces dst_vector_face_access,
+ const DataAccessOnFaces src_vector_face_access) const
+{
+ internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
+ InVector,
+ OutVector,
+ CLASS,
+ true>
+ worker(*this,
+ src,
+ dst,
+ false,
+ *owning_class,
+ cell_operation,
+ face_operation,
+ boundary_operation,
+ src_vector_face_access,
+ dst_vector_face_access,
+ operation_before_loop,
+ operation_after_loop,
+ dof_handler_index_pre_post);
+ task_info.loop(worker);
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+template <typename CLASS, typename OutVector, typename InVector>
+inline void
+MatrixFree<dim, Number, VectorizedArrayType>::loop(
+ void (CLASS::*cell_operation)(const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &),
+ void (CLASS::*face_operation)(const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &),
+ void (CLASS::*boundary_operation)(
+ const MatrixFree &,
+ OutVector &,
+ const InVector &,
+ const std::pair<unsigned int, unsigned int> &),
+ const CLASS * owning_class,
+ OutVector & dst,
+ const InVector &src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_loop,
+ const std::function<void(const unsigned int, const unsigned int)>
+ & operation_after_loop,
+ const unsigned int dof_handler_index_pre_post,
+ const DataAccessOnFaces dst_vector_face_access,
+ const DataAccessOnFaces src_vector_face_access) const
+{
+ internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
+ InVector,
+ OutVector,
+ CLASS,
+ false>
+ worker(*this,
+ src,
+ dst,
+ false,
+ *owning_class,
+ cell_operation,
+ face_operation,
+ boundary_operation,
+ src_vector_face_access,
+ dst_vector_face_access,
+ operation_before_loop,
+ operation_after_loop,
+ dof_handler_index_pre_post);
+ task_info.loop(worker);
+}
+
+
+
template <int dim, typename Number, typename VectorizedArrayType>
template <typename CLASS, typename OutVector, typename InVector>
inline void
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 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 the correctness of the matrix-free loop with additional
+// operation_before_loop and operation_after_loop lambdas
+
+// similar to pre_and_post_loop_01 but uses loop instead of cell_loop with
+// do-nothing face and boundary operations
+
+#include <deal.II/distributed/tria.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/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ typename Number = double>
+class Matrix
+{
+public:
+ Matrix(const MatrixFree<dim, Number> &data_in)
+ : data(data_in)
+ {}
+
+ const std::function<void(const MatrixFree<dim, Number> &,
+ LinearAlgebra::distributed::Vector<Number> &,
+ const LinearAlgebra::distributed::Vector<Number> &,
+ const std::pair<unsigned int, unsigned int> &)>
+ do_nothing = [](const MatrixFree<dim, Number> &,
+ LinearAlgebra::distributed::Vector<Number> &,
+ const LinearAlgebra::distributed::Vector<Number> &,
+ const std::pair<unsigned int, unsigned int> &) {};
+
+ void
+ vmult(LinearAlgebra::distributed::Vector<Number> & dst,
+ const LinearAlgebra::distributed::Vector<Number> &src) const
+ {
+ const std::function<void(const MatrixFree<dim, Number> &,
+ LinearAlgebra::distributed::Vector<Number> &,
+ const LinearAlgebra::distributed::Vector<Number> &,
+ const std::pair<unsigned int, unsigned int> &)>
+ wrap = helmholtz_operator<dim,
+ fe_degree,
+ LinearAlgebra::distributed::Vector<Number>,
+ n_q_points_1d>;
+
+ dst = 0;
+ data.loop(wrap, do_nothing, do_nothing, dst, src);
+ for (auto i : data.get_constrained_dofs())
+ dst.local_element(i) = src.local_element(i);
+ }
+
+ void
+ vmult_with_update_basic(
+ LinearAlgebra::distributed::Vector<Number> &dst,
+ LinearAlgebra::distributed::Vector<Number> &src,
+ LinearAlgebra::distributed::Vector<Number> &other_vector) const
+ {
+ src.add(1.5, other_vector);
+ other_vector.add(0.7, dst);
+ dst = 0;
+ vmult(dst, src);
+ dst.sadd(0.63, -1.3, other_vector);
+ }
+
+ void
+ vmult_with_update_merged(
+ LinearAlgebra::distributed::Vector<Number> &dst,
+ LinearAlgebra::distributed::Vector<Number> &src,
+ LinearAlgebra::distributed::Vector<Number> &other_vector) const
+ {
+ const std::function<void(const MatrixFree<dim, Number> &,
+ LinearAlgebra::distributed::Vector<Number> &,
+ const LinearAlgebra::distributed::Vector<Number> &,
+ const std::pair<unsigned int, unsigned int> &)>
+ wrap = helmholtz_operator<dim,
+ fe_degree,
+ LinearAlgebra::distributed::Vector<Number>,
+ n_q_points_1d>;
+ data.loop(
+ wrap,
+ do_nothing,
+ do_nothing,
+ dst,
+ src,
+ // operation before cell operation
+ [&](const unsigned int start_range, const unsigned int end_range) {
+ for (unsigned int i = start_range; i < end_range; ++i)
+ {
+ src.local_element(i) += 1.5 * other_vector.local_element(i);
+ other_vector.local_element(i) += 0.7 * dst.local_element(i);
+ dst.local_element(i) = 0;
+ }
+ },
+ // operation after cell operation
+ [&](const unsigned int start_range, const unsigned int end_range) {
+ for (unsigned int i = start_range; i < end_range; ++i)
+ {
+ dst.local_element(i) =
+ 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i);
+ }
+ });
+ }
+
+private:
+ const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(7 - dim);
+
+ FE_Q<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ MatrixFree<dim, number> mf_data;
+ {
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+ data.mapping_update_flags_inner_faces = update_JxW_values;
+ mf_data.reinit(MappingQ1<dim>{}, dof, constraints, quad, data);
+ }
+
+ Matrix<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+ LinearAlgebra::distributed::Vector<number> vec1, vec2, vec3;
+ mf_data.initialize_dof_vector(vec1);
+ mf_data.initialize_dof_vector(vec2);
+ mf_data.initialize_dof_vector(vec3);
+
+ for (unsigned int i = 0; i < vec1.local_size(); ++i)
+ {
+ // Multiply by 0.01 to make float error with roundoff less than the
+ // numdiff absolute tolerance
+ double entry = 0.01 * random_value<double>();
+ vec1.local_element(i) = entry;
+ entry = 0.01 * random_value<double>();
+ vec2.local_element(i) = entry;
+ entry = 0.01 * random_value<double>();
+ vec3.local_element(i) = entry;
+ }
+
+ LinearAlgebra::distributed::Vector<number> ref1 = vec1;
+ LinearAlgebra::distributed::Vector<number> ref2 = vec2;
+ LinearAlgebra::distributed::Vector<number> ref3 = vec3;
+
+ mf.vmult_with_update_basic(ref3, ref2, ref1);
+ mf.vmult_with_update_merged(vec3, vec2, vec1);
+
+ vec3 -= ref3;
+ deallog << "Error in 1x merged operation: " << vec3.linfty_norm()
+ << std::endl;
+
+ ref3 = 0;
+
+ mf.vmult_with_update_basic(ref1, ref2, ref3);
+ mf.vmult_with_update_merged(vec1, vec2, vec3);
+
+ mf.vmult_with_update_basic(ref1, ref2, ref3);
+ mf.vmult_with_update_merged(vec1, vec2, vec3);
+
+ vec3 -= ref3;
+ deallog << "Error in 2x merged operation: " << vec3.linfty_norm()
+ << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+ mpi_initlog();
+ deallog << std::setprecision(3);
+
+ test<2, 3, double>();
+ test<2, 3, float>();
+ test<3, 3, double>();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// perform some boundary and face operation but fully overwrite the result
+// within the post operation
+
+#include <deal.II/distributed/tria.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/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ typename Number = double>
+class Matrix
+{
+public:
+ Matrix(const MatrixFree<dim, Number> &data_in)
+ : data(data_in)
+ {}
+
+ void
+ do_nothing(const MatrixFree<dim, Number> &,
+ LinearAlgebra::distributed::Vector<Number> &,
+ const LinearAlgebra::distributed::Vector<Number> &,
+ const std::pair<unsigned int, unsigned int> &) const {};
+
+ void
+ do_boundary(const MatrixFree<dim, Number> & mf,
+ LinearAlgebra::distributed::Vector<Number> & dst,
+ const LinearAlgebra::distributed::Vector<Number> &src,
+ const std::pair<unsigned int, unsigned int> & range) const
+ {
+ using FEFaceIntegrator = FEFaceEvaluation<dim, -1, 0, 1, Number>;
+ FEFaceIntegrator eval(mf, true);
+
+ for (unsigned int face = range.first; face < range.second; ++face)
+ {
+ eval.reinit(face);
+ eval.gather_evaluate(src, EvaluationFlags::values);
+ for (unsigned int q = 0; q < eval.n_q_points; ++q)
+ {
+ eval.submit_value(1.23 * eval.get_value(q), q);
+ }
+ eval.integrate_scatter(EvaluationFlags::values, dst);
+ }
+ };
+
+ void
+ do_faces(const MatrixFree<dim, Number> & mf,
+ LinearAlgebra::distributed::Vector<Number> & dst,
+ const LinearAlgebra::distributed::Vector<Number> &src,
+ const std::pair<unsigned int, unsigned int> & range) const
+ {
+ using FEFaceIntegrator = FEFaceEvaluation<dim, -1, 0, 1, Number>;
+ FEFaceIntegrator eval(mf, true);
+
+ for (unsigned int face = range.first; face < range.second; ++face)
+ {
+ eval.reinit(face);
+ eval.gather_evaluate(src, EvaluationFlags::values);
+ for (unsigned int q = 0; q < eval.n_q_points; ++q)
+ {
+ eval.submit_value(1.23 * eval.get_value(q), q);
+ }
+ eval.integrate_scatter(EvaluationFlags::values, dst);
+ }
+ };
+
+ void
+ vmult_pre_post_op_only(
+ LinearAlgebra::distributed::Vector<Number> & dst,
+ const LinearAlgebra::distributed::Vector<Number> &src) const
+ {
+ data.loop(
+ &Matrix::do_nothing,
+ &Matrix::do_faces,
+ &Matrix::do_boundary,
+ this,
+ dst,
+ src,
+ // operation before
+ [&](const unsigned int start_range, const unsigned int end_range) {
+ for (unsigned int i = start_range; i < end_range; ++i)
+ {
+ dst.local_element(i) = 0;
+ }
+ },
+ // operation after
+ [&](const unsigned int start_range, const unsigned int end_range) {
+ for (unsigned int i = start_range; i < end_range; ++i)
+ {
+ dst.local_element(i) = 2 * src.local_element(i);
+ }
+ });
+ }
+
+private:
+ const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(7 - dim);
+
+ FE_Q<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ MatrixFree<dim, number> mf_data;
+ {
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+ data.mapping_update_flags_inner_faces = update_JxW_values;
+ mf_data.reinit(MappingQ1<dim>{}, dof, constraints, quad, data);
+ }
+
+ Matrix<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+ LinearAlgebra::distributed::Vector<number> vec1, vec2, vec3;
+ mf_data.initialize_dof_vector(vec1);
+ mf_data.initialize_dof_vector(vec2);
+
+ for (unsigned int i = 0; i < vec1.local_size(); ++i)
+ {
+ // Multiply by 0.01 to make float error with roundoff less than the
+ // numdiff absolute tolerance
+ double entry = 0.01 * random_value<double>();
+ vec1.local_element(i) = entry;
+ entry = 0.01 * random_value<double>();
+ vec2.local_element(i) = entry;
+ }
+
+ mf.vmult_pre_post_op_only(vec2, vec1);
+
+ auto ref = vec1;
+ ref *= 2.0;
+ vec2 -= ref;
+ deallog << "Error: " << vec2.linfty_norm() << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+ mpi_initlog();
+ deallog << std::setprecision(3);
+
+ test<2, 3, double>();
+ test<2, 3, float>();
+ test<3, 3, double>();
+}
projects / dealii.git / commitdiff