// Simple case: We have contiguous storage, so we can simply copy out the
// data
- if (this->dof_info->index_storage_variants[ind][this->cell] ==
- internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- interleaved_contiguous &&
- n_lanes == VectorizedArrayType::size())
+ if ((this->dof_info->index_storage_variants[ind][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous &&
+ n_lanes == VectorizedArrayType::size()) &&
+ !(is_face &&
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false))
{
const unsigned int dof_index =
dof_indices_cont[this->cell * VectorizedArrayType::size()] +
return;
}
+ std::array<unsigned int, VectorizedArrayType::size()> cells =
+ this->get_cell_or_face_ids();
+
// More general case: Must go through the components one by one and apply
// some transformations
const unsigned int n_filled_lanes =
this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
- unsigned int dof_indices[VectorizedArrayType::size()];
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- dof_indices[v] =
- dof_indices_cont[this->cell * VectorizedArrayType::size() + v] +
- this->dof_info
- ->component_dof_indices_offset[this->active_fe_index]
- [this->first_selected_component] *
- this->dof_info->dof_indices_interleave_strides
- [ind][this->cell * VectorizedArrayType::size() + v];
-
- for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
- dof_indices[v] = numbers::invalid_unsigned_int;
+ const bool is_ecl =
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false;
if (vectors_sm[0] != nullptr)
{
for (unsigned int v = 0; v < n_filled_lanes; ++v)
{
+ Assert(cells[v] != numbers::invalid_unsigned_int,
+ ExcNotImplemented());
Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
ExcIndexRange(
ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
- Assert(this->cell * VectorizedArrayType::size() + v <
+ Assert(cells[v] <
this->dof_info->dof_indices_contiguous_sm[ind].size(),
ExcIndexRange(
- this->cell * VectorizedArrayType::size() + v,
+ cells[v],
0,
this->dof_info->dof_indices_contiguous_sm[ind].size()));
const auto &temp =
- this->dof_info->dof_indices_contiguous_sm
- [ind][this->cell * VectorizedArrayType::size() + v];
+ this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
if (temp.first != numbers::invalid_unsigned_int)
vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
};
if (n_filled_lanes == VectorizedArrayType::size() &&
- n_lanes == VectorizedArrayType::size())
+ n_lanes == VectorizedArrayType::size() && !is_ecl)
{
if (n_components == 1 || n_fe_components == 1)
{
return;
}
+ unsigned int dof_indices[VectorizedArrayType::size()];
+
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ {
+ Assert(cells[v] != numbers::invalid_unsigned_int, ExcNotImplemented());
+ dof_indices[v] =
+ dof_indices_cont[cells[v]] +
+ this->dof_info
+ ->component_dof_indices_offset[this->active_fe_index]
+ [this->first_selected_component] *
+ this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
+ }
+
+ for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
+ dof_indices[v] = numbers::invalid_unsigned_int;
+
// In the case with contiguous cell indices, we know that there are no
// constraints and that the indices within each element are contiguous
if (n_filled_lanes == VectorizedArrayType::size() &&
- n_lanes == VectorizedArrayType::size())
+ n_lanes == VectorizedArrayType::size() && !is_ecl)
{
if (this->dof_info->index_storage_variants[ind][this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
!(evaluation_flag & EvaluationFlags::gradients))
return;
- if (fe_degree > -1)
- internal::FEFaceEvaluationImplEvaluateSelector<dim, VectorizedArrayType>::
- template run<fe_degree, n_q_points_1d>(
- n_components,
- *this->data,
- values_array,
- this->begin_values(),
- this->begin_gradients(),
- this->scratch_data,
- evaluation_flag & EvaluationFlags::values,
- evaluation_flag & EvaluationFlags::gradients,
- this->face_no,
- this->subface_index,
- this->face_orientation,
- this->descriptor->face_orientations);
+ if (this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false)
+ {
+ const auto face_nos = this->compute_face_no_data();
+ const auto face_orientations = this->compute_face_orientations();
+
+# ifdef DEBUG
+ // currently on structured meshes are supported -> face numers and
+ // orientations have to be the same for all filled lanes
+ for (unsigned int v = 1; v < VectorizedArrayType::size(); ++v)
+ {
+ if (face_nos[v] != numbers::invalid_unsigned_int)
+ AssertDimension(face_nos[0], face_nos[v]);
+ if (face_orientations[v] != numbers::invalid_unsigned_int)
+ AssertDimension(face_orientations[0], face_orientations[v]);
+ }
+# endif
+
+ internal::FEFaceEvaluationImplEvaluateSelector<dim, VectorizedArrayType>::
+ template run<fe_degree, n_q_points_1d>(
+ n_components,
+ *this->data,
+ values_array,
+ this->begin_values(),
+ this->begin_gradients(),
+ this->scratch_data,
+ evaluation_flag & EvaluationFlags::values,
+ evaluation_flag & EvaluationFlags::gradients,
+ face_nos[0],
+ this->subface_index,
+ face_orientations[0],
+ this->descriptor->face_orientations);
+ }
else
- internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
- evaluate(n_components,
- *this->data,
- values_array,
- this->begin_values(),
- this->begin_gradients(),
- this->scratch_data,
- evaluation_flag & EvaluationFlags::values,
- evaluation_flag & EvaluationFlags::gradients,
- this->face_no,
- this->subface_index,
- this->face_orientation,
- this->descriptor->face_orientations);
+ {
+ if (fe_degree > -1)
+ internal::FEFaceEvaluationImplEvaluateSelector<dim,
+ VectorizedArrayType>::
+ template run<fe_degree, n_q_points_1d>(
+ n_components,
+ *this->data,
+ values_array,
+ this->begin_values(),
+ this->begin_gradients(),
+ this->scratch_data,
+ evaluation_flag & EvaluationFlags::values,
+ evaluation_flag & EvaluationFlags::gradients,
+ this->face_no,
+ this->subface_index,
+ this->face_orientation,
+ this->descriptor->face_orientations);
+ else
+ internal::FEFaceEvaluationFactory<dim, Number, VectorizedArrayType>::
+ evaluate(n_components,
+ *this->data,
+ values_array,
+ this->begin_values(),
+ this->begin_gradients(),
+ this->scratch_data,
+ evaluation_flag & EvaluationFlags::values,
+ evaluation_flag & EvaluationFlags::gradients,
+ this->face_no,
+ this->subface_index,
+ this->face_orientation,
+ this->descriptor->face_orientations);
+ }
# ifdef DEBUG
if (evaluation_flag & EvaluationFlags::values)
}
else
{
- std::fill(face_no_data.begin(), face_no_data.end(), 0);
+ std::fill(face_no_data.begin(),
+ face_no_data.end(),
+ numbers::invalid_unsigned_int);
if (dim == 3)
{
face_no_data[i] = face_orientation;
}
}
+ else
+ {
+ std::fill(
+ face_no_data.begin(),
+ face_no_data.begin() +
+ this->dof_info->n_vectorization_lanes_filled
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell],
+ 0);
+ }
}
return face_no_data;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+// Compare the evaluation of a SIP Laplace operator with face- and element-
+// centric loops on a structured grid.
+
+template <int dim, typename Number>
+class CosineFunction : public Function<dim, Number>
+{
+public:
+ Number
+ value(const Point<dim> &p, const unsigned int component = 0) const
+ {
+ (void)component;
+ return std::cos(2 * numbers::PI * p[0]);
+ }
+};
+
+template <int dim,
+ int fe_degree,
+ int n_points = fe_degree + 1,
+ typename Number = double,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test(const unsigned int n_refinements = 1)
+{
+ using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria, 0.0, 1.0, true);
+
+ if (false)
+ {
+ std::vector<dealii::GridTools::PeriodicFacePair<
+ typename dealii::Triangulation<dim>::cell_iterator>>
+ periodic_faces;
+
+ if (dim >= 1)
+ dealii::GridTools::collect_periodic_faces(
+ tria, 0, 1, 0, periodic_faces);
+
+ if (dim >= 2)
+ dealii::GridTools::collect_periodic_faces(
+ tria, 2, 3, 1, periodic_faces);
+
+ if (dim >= 3)
+ dealii::GridTools::collect_periodic_faces(
+ tria, 4, 5, 2, periodic_faces);
+
+ tria.add_periodicity(periodic_faces);
+ }
+
+ tria.refine_global(n_refinements);
+
+ FE_DGQ<dim> fe(fe_degree);
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ MappingQ<dim> mapping(1);
+ QGauss<1> quad(n_points);
+
+ AffineConstraints<Number> constraint;
+
+ using MF = MatrixFree<dim, Number, VectorizedArrayType>;
+
+ typename MF::AdditionalData additional_data;
+ additional_data.mapping_update_flags = update_values | update_gradients;
+ additional_data.mapping_update_flags_inner_faces =
+ update_values | update_gradients;
+ additional_data.mapping_update_flags_boundary_faces =
+ update_values | update_gradients;
+ additional_data.tasks_parallel_scheme =
+ MF::AdditionalData::TasksParallelScheme::none;
+ // additional_data.mapping_update_flags_faces_by_cells = update_values;
+ // additional_data.hold_all_faces_to_owned_cells = true;
+
+ MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data);
+
+ MF matrix_free;
+ matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+ VectorType src, dst;
+
+ matrix_free.initialize_dof_vector(src);
+ matrix_free.initialize_dof_vector(dst);
+
+ FEEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType> phi(
+ matrix_free);
+ FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+ phi_m(matrix_free, true);
+ FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+ phi_p(matrix_free, false);
+
+
+ CosineFunction<dim, Number> function;
+ VectorTools::interpolate(dof_handler, function, src);
+
+ dst = 0.0;
+
+ /**
+ * Face-centric loop
+ */
+ matrix_free.template loop<VectorType, VectorType>(
+ [&](const auto &, auto &dst, const auto &src, const auto range) {
+ for (unsigned int cell = range.first; cell < range.second; ++cell)
+ {
+ phi.reinit(cell);
+ phi.read_dof_values(src);
+ phi.evaluate(false, true, false);
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ phi.submit_gradient(phi.get_gradient(q), q);
+ phi.integrate(false, true);
+ phi.set_dof_values(dst);
+ }
+ },
+ [&](const auto &, auto &dst, const auto &src, const auto range) {
+ for (unsigned int face = range.first; face < range.second; face++)
+ {
+ phi_m.reinit(face);
+ phi_p.reinit(face);
+
+ phi_m.read_dof_values(src);
+ phi_m.evaluate(true, true);
+ phi_p.read_dof_values(src);
+ phi_p.evaluate(true, true);
+ VectorizedArrayType sigmaF =
+ (std::abs((phi_m.get_normal_vector(0) *
+ phi_m.inverse_jacobian(0))[dim - 1]) +
+ std::abs((phi_m.get_normal_vector(0) *
+ phi_p.inverse_jacobian(0))[dim - 1])) *
+ (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+ for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+ {
+ VectorizedArrayType average_value =
+ (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
+ VectorizedArrayType average_valgrad =
+ phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q);
+ average_valgrad =
+ average_value * 2. * sigmaF - average_valgrad * 0.5;
+ phi_m.submit_normal_derivative(-average_value, q);
+ phi_p.submit_normal_derivative(-average_value, q);
+ phi_m.submit_value(average_valgrad, q);
+ phi_p.submit_value(-average_valgrad, q);
+ }
+ phi_m.integrate(true, true);
+ phi_m.distribute_local_to_global(dst);
+ phi_p.integrate(true, true);
+ phi_p.distribute_local_to_global(dst);
+ }
+ },
+ [&](const auto &, auto &dst, const auto &src, const auto face_range) {
+ for (unsigned int face = face_range.first; face < face_range.second;
+ face++)
+ {
+ phi_m.reinit(face);
+ phi_m.read_dof_values(src);
+ phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+ VectorizedArrayType sigmaF =
+ std::abs((phi_m.get_normal_vector(0) *
+ phi_m.inverse_jacobian(0))[dim - 1]) *
+ Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+ for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+ {
+ VectorizedArrayType average_value = phi_m.get_value(q);
+ VectorizedArrayType average_valgrad =
+ -phi_m.get_normal_derivative(q);
+ average_valgrad += average_value * sigmaF * 2.0;
+ phi_m.submit_normal_derivative(-average_value, q);
+ phi_m.submit_value(average_valgrad, q);
+ }
+
+ phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+ phi_m.distribute_local_to_global(dst);
+ }
+ },
+ dst,
+ src,
+ false,
+ MF::DataAccessOnFaces::gradients,
+ MF::DataAccessOnFaces::gradients);
+
+ deallog << dst.l2_norm() << std::endl;
+ dst.print(deallog.get_file_stream());
+
+ dst = 0.0;
+
+ /**
+ * Element-centric loop
+ */
+ matrix_free.template loop_cell_centric<VectorType, VectorType>(
+ [&](const auto &, auto &dst, const auto &src, const auto range) {
+ for (unsigned int cell = range.first; cell < range.second; ++cell)
+ {
+ phi.reinit(cell);
+ phi.read_dof_values(src);
+ phi.evaluate(false, true, false);
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ phi.submit_gradient(phi.get_gradient(q), q);
+ phi.integrate(false, true);
+
+ for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ {
+ auto bids =
+ matrix_free.get_faces_by_cells_boundary_id(cell, face);
+
+ if (bids[0] != numbers::internal_face_boundary_id)
+ {
+ phi_m.reinit(cell, face);
+ phi_m.read_dof_values(src);
+ phi_m.evaluate(EvaluationFlags::values |
+ EvaluationFlags::gradients);
+ VectorizedArrayType sigmaF =
+ std::abs((phi_m.get_normal_vector(0) *
+ phi_m.inverse_jacobian(0))[dim - 1]) *
+ Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+ for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+ {
+ VectorizedArrayType average_value = phi_m.get_value(q);
+ VectorizedArrayType average_valgrad =
+ -phi_m.get_normal_derivative(q);
+ average_valgrad += average_value * sigmaF * 2.0;
+ phi_m.submit_normal_derivative(-average_value, q);
+ phi_m.submit_value(average_valgrad, q);
+ }
+
+ phi_m.integrate(EvaluationFlags::values |
+ EvaluationFlags::gradients);
+
+ for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+ phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q];
+ }
+ else
+ {
+ phi_m.reinit(cell, face);
+ phi_p.reinit(cell, face);
+
+ phi_m.read_dof_values(src);
+ phi_m.evaluate(EvaluationFlags::values |
+ EvaluationFlags::gradients);
+ phi_p.read_dof_values(src);
+ phi_p.evaluate(EvaluationFlags::values |
+ EvaluationFlags::gradients);
+
+ VectorizedArrayType sigmaF =
+ (std::abs((phi_m.get_normal_vector(0) *
+ phi_m.inverse_jacobian(0))[dim - 1]) +
+ std::abs((phi_m.get_normal_vector(0) *
+ phi_p.inverse_jacobian(0))[dim - 1])) *
+ (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+ for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+ {
+ VectorizedArrayType average_value =
+ (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
+ VectorizedArrayType average_valgrad =
+ phi_m.get_normal_derivative(q) +
+ phi_p.get_normal_derivative(q);
+ average_valgrad =
+ average_value * 2. * sigmaF - average_valgrad * 0.5;
+ phi_m.submit_normal_derivative(-average_value, q);
+ phi_m.submit_value(average_valgrad, q);
+ }
+ phi_m.integrate(EvaluationFlags::values |
+ EvaluationFlags::gradients);
+
+ for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+ phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q];
+ }
+ }
+
+ phi.set_dof_values(dst);
+ }
+ },
+ dst,
+ src,
+ false,
+ MF::DataAccessOnFaces::gradients);
+
+ deallog << dst.l2_norm() << std::endl;
+ dst.print(deallog.get_file_stream());
+}
+
+int
+main()
+{
+ initlog();
+ test<2, 2, 3, double, VectorizedArray<double>>();
+}
projects / dealii.git / commitdiff