#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/matrix_free/type_traits.h>
DEAL_II_NAMESPACE_OPEN
data, temp1, values_array, integrate_gradients, face_no);
}
- template <std::size_t n_face_orientations>
+ template <std::size_t n_face_orientations,
+ typename VectorType,
+ typename std::enable_if<
+ !internal::has_begin<VectorType>::value ||
+ !std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number2 *>::value,
+ VectorType>::type * = nullptr>
+ static bool
+ gather_evaluate(const VectorType &,
+ const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &,
+ const MatrixFreeFunctions::DoFInfo &,
+ VectorizedArrayType *,
+ VectorizedArrayType *,
+ VectorizedArrayType *,
+ const bool,
+ const bool,
+ const unsigned int,
+ const unsigned int,
+ const std::array<unsigned int, n_face_orientations>,
+ const std::array<unsigned int, n_face_orientations>,
+ const unsigned int,
+ const MatrixFreeFunctions::DoFInfo::DoFAccessIndex,
+ const std::array<unsigned int, n_face_orientations>,
+ const Table<2, unsigned int> &)
+ {
+ return false;
+ }
+
+ template <std::size_t n_face_orientations,
+ typename VectorType,
+ typename std::enable_if<
+ internal::has_begin<VectorType>::value &&
+ std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number2 *>::value,
+ VectorType>::type * = nullptr>
static bool
gather_evaluate(
- const Number2 * src_ptr,
+ const VectorType & src,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
const MatrixFreeFunctions::DoFInfo & dof_info,
VectorizedArrayType * values_quad,
const std::array<unsigned int, n_face_orientations> face_orientations,
const Table<2, unsigned int> & orientation_map)
{
+ static_assert((std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number2 *>::value),
+ "Types do not match.");
+
return process_and_io( //
false /*=evaluate*/,
- src_ptr,
+ src.begin(),
data,
dof_info,
values_quad,
});
}
- template <std::size_t n_face_orientations>
+ template <std::size_t n_face_orientations,
+ typename VectorType,
+ typename std::enable_if<
+ !internal::has_begin<VectorType>::value ||
+ !std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number2 *>::value,
+ VectorType>::type * = nullptr>
+ static bool
+ integrate_scatter(
+ VectorType &,
+ const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
+ const MatrixFreeFunctions::DoFInfo &,
+ VectorizedArrayType *values_array,
+ VectorizedArrayType *values_quad,
+ VectorizedArrayType *gradients_quad,
+ VectorizedArrayType *scratch_data,
+ const bool integrate_values,
+ const bool integrate_gradients,
+ const unsigned int,
+ const unsigned int,
+ const std::array<unsigned int, n_face_orientations>,
+ const std::array<unsigned int, n_face_orientations> face_nos,
+ const unsigned int subface_index,
+ const MatrixFreeFunctions::DoFInfo::DoFAccessIndex,
+ const std::array<unsigned int, n_face_orientations> face_orientations,
+ const Table<2, unsigned int> & orientation_map)
+ {
+ AssertDimension(face_nos.size(), 1);
+ AssertDimension(face_orientations.size(), 1);
+
+ // for block vectors simply integrate
+ integrate(data,
+ values_array,
+ values_quad,
+ gradients_quad,
+ scratch_data,
+ integrate_values,
+ integrate_gradients,
+ face_nos[0],
+ subface_index,
+ face_orientations[0],
+ orientation_map);
+
+ // default vector access
+ return false;
+ }
+
+ template <std::size_t n_face_orientations,
+ typename VectorType,
+ typename std::enable_if<
+ internal::has_begin<VectorType>::value &&
+ std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number2 *>::value,
+ VectorType>::type * = nullptr>
static bool
integrate_scatter(
- Number2 * dst_ptr,
+ VectorType & dst,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
const MatrixFreeFunctions::DoFInfo & dof_info,
VectorizedArrayType * values_array,
const std::array<unsigned int, n_face_orientations> face_orientations,
const Table<2, unsigned int> & orientation_map)
{
+ static_assert((std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number2 *>::value),
+ "Types do not match.");
+
return process_and_io( //
true /*=integrate*/,
- dst_ptr,
+ dst.begin(),
data,
dof_info,
values_quad,
gather_evaluate(input_vector, flag);
}
+
+namespace internal
+{
+ /**
+ * Implementation for vectors that have the begin() methods.
+ */
+ template <typename Number,
+ typename VectorizedArrayType,
+ typename VectorType,
+ typename T,
+ typename std::enable_if<
+ internal::has_begin<VectorType>::value &&
+ std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number *>::value,
+ VectorType>::type * = nullptr>
+ bool
+ try_gather_evaluate_inplace(
+ T phi,
+ const VectorType & input_vector,
+ const unsigned int cell,
+ const unsigned int active_fe_index,
+ const unsigned int first_selected_component,
+ const internal::MatrixFreeFunctions::DoFInfo *dof_info,
+ const EvaluationFlags::EvaluationFlags evaluation_flag)
+ {
+ // If the index storage is interleaved and contiguous and the vector storage
+ // has the correct alignment, we can directly pass the pointer into the
+ // vector to the evaluate() call, without reading the vector entries into a
+ // separate data field. This saves some operations.
+ if (std::is_same<typename VectorType::value_type, Number>::value &&
+ dof_info->index_storage_variants
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous &&
+ reinterpret_cast<std::size_t>(
+ input_vector.begin() +
+ dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [cell * VectorizedArrayType::size()]) %
+ sizeof(VectorizedArrayType) ==
+ 0)
+ {
+ const VectorizedArrayType *vec_values =
+ reinterpret_cast<const VectorizedArrayType *>(
+ input_vector.begin() +
+ dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [cell * VectorizedArrayType::size()] +
+ dof_info->component_dof_indices_offset[active_fe_index]
+ [first_selected_component] *
+ VectorizedArrayType::size());
+
+ phi->evaluate(vec_values, evaluation_flag);
+
+ return true;
+ }
+
+ return false;
+ }
+
+ /**
+ * Implementation for all other vectors like block vectors.
+ */
+ template <typename Number,
+ typename VectorizedArrayType,
+ typename VectorType,
+ typename T,
+ typename std::enable_if<
+ !internal::has_begin<VectorType>::value ||
+ !std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number *>::value,
+ VectorType>::type * = nullptr>
+ bool
+ try_gather_evaluate_inplace(T,
+ const VectorType &,
+ const unsigned int,
+ const unsigned int,
+ const unsigned int,
+ const internal::MatrixFreeFunctions::DoFInfo *,
+ const EvaluationFlags::EvaluationFlags)
+ {
+ return false;
+ }
+
+ /**
+ * Implementation for vectors that have the begin() methods.
+ */
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ typename Number,
+ typename VectorizedArrayType,
+ typename VectorType,
+ typename std::enable_if<
+ internal::has_begin<VectorType>::value &&
+ std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number *>::value,
+ VectorType>::type * = nullptr>
+ bool
+ try_integrate_scatter_inplace(
+ VectorType & destination,
+ const unsigned int cell,
+ const unsigned int n_components,
+ const unsigned int active_fe_index,
+ const unsigned int first_selected_component,
+ const internal::MatrixFreeFunctions::DoFInfo *dof_info,
+ VectorizedArrayType * values_quad,
+ VectorizedArrayType * gradients_quad,
+ VectorizedArrayType * scratch_data,
+ const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *data,
+ const EvaluationFlags::EvaluationFlags integration_flag)
+ {
+ // If the index storage is interleaved and contiguous and the vector storage
+ // has the correct alignment, we can directly pass the pointer into the
+ // vector to the integrate() call, without writing temporary results into a
+ // separate data field that will later be added into the vector. This saves
+ // some operations.
+ if (std::is_same<typename VectorType::value_type, Number>::value &&
+ dof_info->index_storage_variants
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous &&
+ reinterpret_cast<std::size_t>(
+ destination.begin() +
+ dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [cell * VectorizedArrayType::size()]) %
+ sizeof(VectorizedArrayType) ==
+ 0)
+ {
+ VectorizedArrayType *vec_values =
+ reinterpret_cast<VectorizedArrayType *>(
+ destination.begin() +
+ dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [cell * VectorizedArrayType::size()] +
+ dof_info->component_dof_indices_offset[active_fe_index]
+ [first_selected_component] *
+ VectorizedArrayType::size());
+ SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
+ integrate(n_components,
+ integration_flag,
+ *data,
+ vec_values,
+ values_quad,
+ gradients_quad,
+ scratch_data,
+ true);
+ return true;
+ }
+
+ return false;
+ }
+
+ /**
+ * Implementation for all other vectors like block vectors.
+ */
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ typename Number,
+ typename VectorizedArrayType,
+ typename VectorType,
+ typename std::enable_if<
+ !internal::has_begin<VectorType>::value ||
+ !std::is_same<decltype(std::declval<VectorType>().begin()),
+ Number *>::value,
+ VectorType>::type * = nullptr>
+ bool
+ try_integrate_scatter_inplace(
+ VectorType &,
+ const unsigned int,
+ const unsigned int,
+ const unsigned int,
+ const unsigned int,
+ const internal::MatrixFreeFunctions::DoFInfo *,
+ const VectorizedArrayType *,
+ const VectorizedArrayType *,
+ const VectorizedArrayType *,
+ const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *,
+ const EvaluationFlags::EvaluationFlags)
+ {
+ return false;
+ }
+} // namespace internal
+
template <int dim,
int fe_degree,
int n_q_points_1d,
gather_evaluate(const VectorType & input_vector,
const EvaluationFlags::EvaluationFlags evaluation_flag)
{
- // If the index storage is interleaved and contiguous and the vector storage
- // has the correct alignment, we can directly pass the pointer into the
- // vector to the evaluate() call, without reading the vector entries into a
- // separate data field. This saves some operations.
- if (std::is_same<typename VectorType::value_type, Number>::value &&
- this->dof_info->index_storage_variants
- [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
- [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
- IndexStorageVariants::interleaved_contiguous &&
- reinterpret_cast<std::size_t>(
- input_vector.begin() +
- this->dof_info->dof_indices_contiguous
- [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
- [this->cell * VectorizedArrayType::size()]) %
- sizeof(VectorizedArrayType) ==
- 0)
- {
- const VectorizedArrayType *vec_values =
- reinterpret_cast<const VectorizedArrayType *>(
- input_vector.begin() +
- this->dof_info->dof_indices_contiguous
- [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
- [this->cell * VectorizedArrayType::size()] +
- this->dof_info
- ->component_dof_indices_offset[this->active_fe_index]
- [this->first_selected_component] *
- VectorizedArrayType::size());
-
- evaluate(vec_values, evaluation_flag);
- }
- else
+ if (internal::try_gather_evaluate_inplace<Number, VectorizedArrayType>(
+ this,
+ input_vector,
+ this->cell,
+ this->active_fe_index,
+ this->first_selected_component,
+ this->dof_info,
+ evaluation_flag) == false)
{
this->read_dof_values(input_vector);
evaluate(this->begin_dof_values(), evaluation_flag);
integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
VectorType & destination)
{
- // If the index storage is interleaved and contiguous and the vector storage
- // has the correct alignment, we can directly pass the pointer into the
- // vector to the integrate() call, without writing temporary results into a
- // separate data field that will later be added into the vector. This saves
- // some operations.
- if (std::is_same<typename VectorType::value_type, Number>::value &&
- this->dof_info->index_storage_variants
- [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
- [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
- IndexStorageVariants::interleaved_contiguous &&
- reinterpret_cast<std::size_t>(
- destination.begin() +
- this->dof_info->dof_indices_contiguous
- [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
- [this->cell * VectorizedArrayType::size()]) %
- sizeof(VectorizedArrayType) ==
- 0)
- {
- VectorizedArrayType *vec_values = reinterpret_cast<VectorizedArrayType *>(
- destination.begin() +
- this->dof_info->dof_indices_contiguous
- [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
- [this->cell * VectorizedArrayType::size()] +
- this->dof_info
- ->component_dof_indices_offset[this->active_fe_index]
- [this->first_selected_component] *
- VectorizedArrayType::size());
- SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
- integrate(n_components,
- integration_flag,
- *this->data,
- vec_values,
- this->values_quad,
- this->gradients_quad,
- this->scratch_data,
- true);
- }
- else
+ if (internal::try_integrate_scatter_inplace<dim,
+ fe_degree,
+ n_q_points_1d,
+ Number,
+ VectorizedArrayType>(
+ destination,
+ this->cell,
+ n_components,
+ this->active_fe_index,
+ this->first_selected_component,
+ this->dof_info,
+ this->values_quad,
+ this->gradients_quad,
+ this->scratch_data,
+ this->data,
+ integration_flag) == false)
{
integrate(integration_flag, this->begin_dof_values());
this->distribute_local_to_global(destination);
gather_evaluate(const VectorType & input_vector,
const EvaluationFlags::EvaluationFlags evaluation_flag)
{
- static_assert(internal::has_begin<VectorType>::value &&
- (std::is_same<decltype(std::declval<VectorType>().begin()),
- double *>::value ||
- std::is_same<decltype(std::declval<VectorType>().begin()),
- float *>::value),
- "This function requires a vector type with begin() function "
- "evaluating to a pointer to basic number (float,double). "
- "Use read_dof_values() followed by evaluate() instead.");
-
Assert(
(evaluation_flag &
~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
n_components,
Number,
VectorizedArrayType>::
- gather_evaluate(input_vector.begin(),
+ gather_evaluate(input_vector,
*this->data,
*this->dof_info,
this->begin_values(),
n_components,
Number,
VectorizedArrayType>::
- gather_evaluate(input_vector.begin(),
+ gather_evaluate(input_vector,
*this->data,
*this->dof_info,
this->begin_values(),
integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
VectorType & destination)
{
- static_assert(internal::has_begin<VectorType>::value &&
- (std::is_same<decltype(std::declval<VectorType>().begin()),
- double *>::value ||
- std::is_same<decltype(std::declval<VectorType>().begin()),
- float *>::value),
- "This function requires a vector type with begin() function "
- "evaluating to a pointer to basic number (float,double). "
- "Use integrate() followed by distribute_local_to_global() "
- "instead.");
-
Assert((this->dof_access_index ==
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
this->is_interior_face == false) == false,
n_components,
Number,
VectorizedArrayType>::
- integrate_scatter(destination.begin(),
+ integrate_scatter(destination,
*this->data,
*this->dof_info,
this->begin_dof_values(),
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// similar to matrix_vector_faces_13 (renumbering of degrees of freedom to
+// better vectorized access, gather_evaluate and integrate_scatter) but with
+// block vectors
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/lac/la_parallel_block_vector.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_faces_common.h"
+
+
+
+template <int dim, int fe_degree_>
+void
+test()
+{
+ // raise element degree by one to test quadratic and cubic shape functions
+ // rather than linears and quadratics according to the main function in
+ // matrix_vector_faces_common.h
+
+ const unsigned int fe_degree = fe_degree_ + 1;
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria, -1, 1);
+
+ for (typename Triangulation<dim>::cell_iterator cell = tria.begin();
+ cell != tria.end();
+ ++cell)
+ for (const unsigned int f : GeometryInfo<dim>::face_indices())
+ if (cell->at_boundary(f))
+ cell->face(f)->set_all_boundary_ids(f);
+ std::vector<
+ GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
+ periodic_faces;
+ for (unsigned int d = 0; d < dim; ++d)
+ GridTools::collect_periodic_faces(
+ tria, 2 * d, 2 * d + 1, d, periodic_faces);
+ tria.add_periodicity(periodic_faces);
+
+ tria.refine_global(3 - dim);
+
+ FE_DGQHermite<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ for (unsigned int test = 0; test < 2; ++test)
+ {
+ tria.refine_global(1);
+
+ dof.distribute_dofs(fe);
+
+ deallog << "Testing " << dof.get_fe().get_name();
+ deallog << " on " << dof.n_dofs() << " DoFs";
+ deallog << std::endl;
+
+ MappingQGeneric<dim> mapping(dof.get_fe().degree + 1);
+
+ LinearAlgebra::distributed::BlockVector<double> in(1);
+ LinearAlgebra::distributed::BlockVector<double> out(1);
+ LinearAlgebra::distributed::BlockVector<double> out_dist(1);
+
+ MatrixFree<dim, double> mf_data;
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, double>::AdditionalData data;
+ data.tasks_parallel_scheme =
+ MatrixFree<dim, double>::AdditionalData::none;
+ data.mapping_update_flags_inner_faces =
+ (update_gradients | update_JxW_values);
+ data.mapping_update_flags_boundary_faces =
+ (update_gradients | update_JxW_values);
+ data.initialize_mapping = false;
+
+ mf_data.reinit(mapping, dof, constraints, quad, data);
+
+ std::vector<types::global_dof_index> renumbering;
+ mf_data.renumber_dofs(renumbering);
+ dof.renumber_dofs(renumbering);
+
+ data.initialize_mapping = true;
+ mf_data.reinit(mapping, dof, constraints, quad, data);
+
+ mf_data.initialize_dof_vector(in.block(0));
+ mf_data.initialize_dof_vector(out.block(0));
+ mf_data.initialize_dof_vector(out_dist.block(0));
+
+ // Set random seed for reproducibility
+ Testing::srand(42);
+ for (unsigned int i = 0; i < in.block(0).local_size(); ++i)
+ {
+ const double entry = Testing::rand() / (double)RAND_MAX;
+ in.block(0).local_element(i) = entry;
+ }
+
+ MatrixFreeTest<dim,
+ fe_degree,
+ fe_degree + 1,
+ double,
+ LinearAlgebra::distributed::BlockVector<double>>
+ mf(mf_data);
+ mf.vmult(out, in);
+
+ MatrixFreeVariant<dim,
+ fe_degree,
+ fe_degree + 1,
+ double,
+ LinearAlgebra::distributed::BlockVector<double>>
+ mf2(mf_data);
+ mf2.vmult(out_dist, in);
+
+ out_dist -= out;
+
+ double diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+ deallog << "Norm of difference: " << diff_norm << " ";
+
+ // test again, now doing matrix-vector product twice
+ mf2.vmult(out_dist, in);
+ mf2.vmult(out_dist, in);
+ out_dist -= out;
+ diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+ deallog << diff_norm << std::endl;
+ }
+}
projects / dealii.git / commitdiff