*this->dof_info);
}
+ const bool accesses_exterior_dofs =
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face() == false;
+
// Case 2: contiguous indices which use reduced storage of indices and can
// use vectorized load/store operations -> go to separate function
if (this->cell != numbers::invalid_unsigned_int)
bool is_contiguous = true;
// check if exterior cells are not contiguous (ECL case)
- if (is_face && !this->interior_face &&
- (this->dof_access_index ==
- internal::MatrixFreeFunctions::DoFInfo::dof_access_cell))
+ if (accesses_exterior_dofs)
{
const std::array<unsigned int, VectorizedArrayType::size()> &cells =
this->get_cell_ids();
internal::is_vectorizable<VectorType, Number>::value>
vector_selector;
- const bool use_vectorized_path = !(masking_is_active || has_hn_constraints);
+ const bool use_vectorized_path =
+ !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
std::array<VectorizedArrayType *, n_components> values_dofs;
Assert(this->cell != numbers::invalid_unsigned_int, ExcNotImplemented());
+ const bool accesses_exterior_dofs =
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face() == false;
+
// Simple case: We have contiguous storage, so we can simply copy out the
// data
- if ((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) &&
- !(!is_face && !this->is_interior_face()))
+ if (dof_info.index_storage_variants[ind][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous &&
+ n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs)
{
const unsigned int dof_index =
dof_indices_cont[this->cell * VectorizedArrayType::size()] +
const unsigned int n_filled_lanes =
dof_info.n_vectorization_lanes_filled[ind][this->cell];
- const bool is_ecl =
- (this->dof_access_index ==
- internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
- this->is_interior_face() == false) ||
- (!is_face && !this->is_interior_face());
-
if (vectors_sm[0] != nullptr)
{
const auto compute_vector_ptrs = [&](const unsigned int comp) {
};
if (n_filled_lanes == VectorizedArrayType::size() &&
- n_lanes == VectorizedArrayType::size() && !is_ecl)
+ n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs)
{
if (n_components == 1 || this->n_fe_components == 1)
{
// 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() && !is_ecl)
+ n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs)
{
if (dof_info.index_storage_variants[ind][this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
{
for (unsigned int i = 0; i < dofs_per_component; ++i)
operation.process_empty(values_dofs[comp][i]);
- if (dof_info.index_storage_variants[ind][this->cell] ==
- internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- contiguous)
- {
- if (n_components == 1 || this->n_fe_components == 1)
- {
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- if (mask[v] == true)
- for (unsigned int i = 0; i < dofs_per_component; ++i)
- operation.process_dof(dof_indices[v] + i,
- *src[comp],
- values_dofs[comp][i][v]);
- }
- else
- {
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- if (mask[v] == true)
- for (unsigned int i = 0; i < dofs_per_component; ++i)
- operation.process_dof(dof_indices[v] + i +
- comp * dofs_per_component,
- *src[0],
- values_dofs[comp][i][v]);
- }
- }
- else
+ if (accesses_exterior_dofs)
{
- const unsigned int *offsets =
- &dof_info.dof_indices_interleave_strides
- [ind][VectorizedArrayType::size() * this->cell];
for (unsigned int v = 0; v < n_filled_lanes; ++v)
- AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
- if (n_components == 1 || this->n_fe_components == 1)
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ if (mask[v] == true)
{
- if (mask[v] == true)
- for (unsigned int i = 0; i < dofs_per_component; ++i)
- operation.process_dof(dof_indices[v] + i * offsets[v],
- *src[comp],
- values_dofs[comp][i][v]);
+ if (dof_info.index_storage_variants
+ [ind][cells[v] / VectorizedArrayType::size()] ==
+ internal::MatrixFreeFunctions::DoFInfo::
+ IndexStorageVariants::contiguous)
+ {
+ if (n_components == 1 || this->n_fe_components == 1)
+ {
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(dof_indices[v] + i,
+ *src[comp],
+ values_dofs[comp][i][v]);
+ }
+ else
+ {
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(dof_indices[v] + i +
+ comp * dofs_per_component,
+ *src[0],
+ values_dofs[comp][i][v]);
+ }
+ }
+ else
+ {
+ const unsigned int offset =
+ dof_info.dof_indices_interleave_strides[ind][cells[v]];
+ AssertIndexRange(offset, VectorizedArrayType::size() + 1);
+ if (n_components == 1 || this->n_fe_components == 1)
+ {
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(dof_indices[v] + i * offset,
+ *src[comp],
+ values_dofs[comp][i][v]);
+ }
+ else
+ {
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(
+ dof_indices[v] +
+ (i + comp * dofs_per_component) * offset,
+ *src[0],
+ values_dofs[comp][i][v]);
+ }
+ }
}
+ }
+ else
+ {
+ if (dof_info.index_storage_variants[ind][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ contiguous)
+ {
+ if (n_components == 1 || this->n_fe_components == 1)
+ {
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ if (mask[v] == true)
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(dof_indices[v] + i,
+ *src[comp],
+ values_dofs[comp][i][v]);
+ }
+ else
+ {
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ if (mask[v] == true)
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(dof_indices[v] + i +
+ comp * dofs_per_component,
+ *src[0],
+ values_dofs[comp][i][v]);
+ }
+ }
else
{
+ const unsigned int *offsets =
+ &dof_info.dof_indices_interleave_strides
+ [ind][VectorizedArrayType::size() * this->cell];
for (unsigned int v = 0; v < n_filled_lanes; ++v)
- if (mask[v] == true)
- for (unsigned int i = 0; i < dofs_per_component; ++i)
- operation.process_dof(dof_indices[v] +
- (i + comp * dofs_per_component) *
- offsets[v],
- *src[0],
- values_dofs[comp][i][v]);
+ AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
+ if (n_components == 1 || this->n_fe_components == 1)
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ {
+ if (mask[v] == true)
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(dof_indices[v] + i * offsets[v],
+ *src[comp],
+ values_dofs[comp][i][v]);
+ }
+ else
+ {
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ if (mask[v] == true)
+ for (unsigned int i = 0; i < dofs_per_component; ++i)
+ operation.process_dof(
+ dof_indices[v] +
+ (i + comp * dofs_per_component) * offsets[v],
+ *src[0],
+ values_dofs[comp][i][v]);
+ }
}
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 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/grid/grid_tools.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"
+
+// Adapted from ecl.h
+
+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 geometry = 0,
+ const unsigned int n_refinements = 1,
+ const bool print_vector = true,
+ const MPI_Comm comm = MPI_COMM_SELF)
+{
+ using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+ if (geometry == 0)
+ GridGenerator::hyper_cube(tria, 0.0, 1.0, true);
+ else if (geometry == 1)
+ GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1.0);
+ else if (geometry == 2)
+ GridGenerator::hyper_ball(tria);
+ else
+ Assert(false, ExcNotImplemented());
+
+ tria.reset_all_manifolds();
+
+ 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_Q<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 | update_gradients;
+ additional_data.hold_all_faces_to_owned_cells = true;
+ additional_data.communicator_sm = comm;
+
+ 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(EvaluationFlags::gradients);
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ phi.submit_gradient(phi.get_gradient(q), q);
+ phi.integrate(EvaluationFlags::gradients);
+ phi.distribute_local_to_global(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(EvaluationFlags::values | EvaluationFlags::gradients);
+ phi_p.read_dof_values(src);
+ phi_p.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+ VectorizedArrayType sigmaF =
+ (std::abs(
+ (phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) +
+ std::abs(
+ (phi_m.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(EvaluationFlags::values | EvaluationFlags::gradients);
+ phi_m.distribute_local_to_global(dst);
+ phi_p.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+ 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.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;
+
+ if (print_vector)
+ dst.print(deallog.get_file_stream(), 10);
+
+ 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(EvaluationFlags::gradients);
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ phi.submit_gradient(phi.get_gradient(q), q);
+ phi.integrate(EvaluationFlags::gradients);
+
+ for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ {
+ const auto bids =
+ matrix_free.get_faces_by_cells_boundary_id(cell, face);
+
+ std::bitset<VectorizedArrayType::size()> mask;
+ for (unsigned int v = 0;
+ v < matrix_free.n_active_entries_per_cell_batch(cell);
+ ++v)
+ if (bids[v] == numbers::internal_face_boundary_id)
+ mask[v] = true;
+
+ // boundary
+ {
+ phi_m.reinit(cell, face);
+ phi_m.read_dof_values(src);
+ phi_m.evaluate(EvaluationFlags::values |
+ EvaluationFlags::gradients);
+ VectorizedArrayType sigmaF =
+ std::abs((phi_m.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 v = 0; v < VectorizedArrayType::size(); ++v)
+ for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+ phi.begin_dof_values()[q][v] +=
+ mask[v] == false ? phi_m.begin_dof_values()[q][v] : 0.;
+ }
+
+ // interior face
+ {
+ 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, 0, mask);
+ phi_p.evaluate(EvaluationFlags::values |
+ EvaluationFlags::gradients);
+
+ VectorizedArrayType sigmaF =
+ (std::abs((phi_m.normal_vector(0) *
+ phi_m.inverse_jacobian(0))[dim - 1]) +
+ std::abs((phi_m.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 v = 0; v < VectorizedArrayType::size(); ++v)
+ for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+ phi.begin_dof_values()[q][v] +=
+ mask[v] == true ? phi_m.begin_dof_values()[q][v] : 0.;
+ }
+ }
+
+ phi.distribute_local_to_global(dst);
+ }
+ },
+ dst,
+ src,
+ false,
+ MF::DataAccessOnFaces::gradients);
+
+ deallog << dst.l2_norm() << std::endl;
+
+ if (print_vector)
+ dst.print(deallog.get_file_stream(), 10);
+}
projects / dealii.git / commitdiff