const AlignedVector<std::array<T, VectorizedArrayType::n_array_elements>>
&array) const;
+ /**
+ * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
+ * associated with.
+ */
+ std::array<unsigned int, VectorizedArrayType::n_array_elements>
+ get_cell_ids() const;
+
//@}
/**
return jacobian[q_index];
}
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::n_array_elements>
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ get_cell_ids() const
+{
+ const unsigned int v_len = VectorizedArrayType::n_array_elements;
+ std::array<unsigned int, VectorizedArrayType::n_array_elements> cells;
+
+ // initialize array
+ for (unsigned int i = 0; i < v_len; ++i)
+ cells[i] = numbers::invalid_unsigned_int;
+
+ if ((is_face == false) ||
+ (is_face &&
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face))
+ {
+ // cell or interior face face (element-centric loop)
+ for (unsigned int i = 0; i < v_len; ++i)
+ cells[i] = cell * v_len + i;
+ }
+ else if (is_face &&
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false)
+ {
+ // exterior face (element-centric loop): for this case, we need to
+ // look into the FaceInfo field that collects information from both
+ // sides of a face once for the global mesh, and pick the face id that
+ // is not the local one (cell_this).
+ for (unsigned int i = 0; i < v_len; i++)
+ {
+ // compute actual (non vectorized) cell ID
+ const unsigned int cell_this = this->cell * v_len + i;
+ // compute face ID
+ unsigned int face_index =
+ this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+ this->face_no,
+ i);
+
+ if (face_index == numbers::invalid_unsigned_int)
+ continue; // invalid face ID: no neighbor on boundary
+
+ // get cell ID on both sides of face
+ auto cell_m = this->matrix_info->get_face_info(face_index / v_len)
+ .cells_interior[face_index % v_len];
+ auto cell_p = this->matrix_info->get_face_info(face_index / v_len)
+ .cells_exterior[face_index % v_len];
+
+ // compare the IDs with the given cell ID
+ if (cell_m == cell_this)
+ cells[i] = cell_p; // neighbor has the other ID
+ else if (cell_p == cell_this)
+ cells[i] = cell_m;
+ }
+ }
+ else if (is_face)
+ {
+ // face-centric faces
+ const unsigned int *cells_ =
+ is_interior_face ?
+ &this->matrix_info->get_face_info(cell).cells_interior[0] :
+ &this->matrix_info->get_face_info(cell).cells_exterior[0];
+ for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i)
+ if (cells_[i] != numbers::invalid_unsigned_int)
+ cells[i] = cells_[i];
+ }
+
+ return cells;
+}
+
template <int dim,
Assert(matrix_info != nullptr, ExcNotImplemented());
AssertDimension(array.size(),
matrix_info->get_task_info().cell_partition_data.back());
- if (is_face)
- {
- VectorizedArrayType out = make_vectorized_array<Number>(Number(1.));
- const unsigned int *cells =
- is_interior_face ?
- &this->matrix_info->get_face_info(cell).cells_interior[0] :
- &this->matrix_info->get_face_info(cell).cells_exterior[0];
- for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i)
- if (cells[i] != numbers::invalid_unsigned_int)
- out[i] = array[cells[i] / VectorizedArrayType::n_array_elements]
- [cells[i] % VectorizedArrayType::n_array_elements];
- return out;
- }
- else
- return array[cell];
+
+ // 1) collect ids of cell
+ const auto cells = this->get_cell_ids();
+
+ // 2) actually gather values
+ VectorizedArrayType out = make_vectorized_array<Number>(Number(1.));
+ for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i)
+ if (cells[i] != numbers::invalid_unsigned_int)
+ out[i] = array[cells[i] / VectorizedArrayType::n_array_elements]
+ [cells[i] % VectorizedArrayType::n_array_elements];
+ return out;
}
Assert(matrix_info != nullptr, ExcNotImplemented());
AssertDimension(array.size(),
matrix_info->get_task_info().cell_partition_data.back());
- if (is_face)
- {
- std::array<T, VectorizedArrayType::n_array_elements> out;
- const unsigned int * cells =
- is_interior_face ?
- &this->matrix_info->get_face_info(cell).cells_interior[0] :
- &this->matrix_info->get_face_info(cell).cells_exterior[0];
- for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i)
- if (cells[i] != numbers::invalid_unsigned_int)
- out[i] = array[cells[i] / VectorizedArrayType::n_array_elements]
- [cells[i] % VectorizedArrayType::n_array_elements];
- return out;
- }
- else
- return array[cell];
+
+ // 1) collect ids of cell
+ const auto cells = this->get_cell_ids();
+
+ // 2) actually gather values
+ std::array<T, VectorizedArrayType::n_array_elements> out;
+ for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i)
+ if (cells[i] != numbers::invalid_unsigned_int)
+ out[i] = array[cells[i] / VectorizedArrayType::n_array_elements]
+ [cells[i] % VectorizedArrayType::n_array_elements];
+ return out;
}
Assert(this->mapped_geometry == nullptr,
ExcMessage("FEEvaluation was initialized without a matrix-free object."
" Integer indexing is not possible"));
- Assert(this->is_interior_face == true,
- ExcMessage(
- "Cell-based FEFaceEvaluation::reinit only possible for the "
- "interior face with second argument to constructor as true"));
if (this->mapped_geometry != nullptr)
return;
Assert(this->matrix_info != nullptr, ExcNotInitialized());
VectorizedArrayType::n_array_elements> &
get_face_info(const unsigned int face_batch_number) const;
+
+ /**
+ * Return the table that translates a triple of the macro cell number,
+ * the index of a face within a cell and the index within the cell batch of
+ * vectorization into the index within the faces array.
+ */
+ const Table<3, unsigned int> &
+ get_cell_and_face_to_plain_faces() const;
+
/**
* Obtains a scratch data object for internal use. Make sure to release it
* afterwards by passing the pointer you obtain from this object to the
+template <int dim, typename Number, typename VectorizedArrayType>
+inline const Table<3, unsigned int> &
+MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
+ const
+{
+ return face_info.cell_and_face_to_plain_faces;
+}
+
+
+
template <int dim, typename Number, typename VectorizedArrayType>
inline const Quadrature<dim> &
MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
true);
face_info.cell_and_face_boundary_id.fill(numbers::invalid_boundary_id);
- for (unsigned int f = 0; f < task_info.boundary_partition_data.back();
+ for (unsigned int f = 0; f < task_info.ghost_face_partition_data.back();
++f)
for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements &&
face_info.faces[f].cells_interior[v] !=
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests matrix-free read_cell_data for each cell (incl. inner and outer faces,
+// which is needed for element-centric-loops)
+
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.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/multigrid/mg_constrained_dofs.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim,
+ typename Number = double,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+class Test
+{
+public:
+ using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+ void
+ run(unsigned int fe_degree)
+ {
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ FE_DGQ<dim> fe(fe_degree);
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ MappingQ<dim> mapping(1);
+
+ QGauss<1> quad(fe_degree + 1);
+
+ AffineConstraints<Number> constraint;
+
+ typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+ additional_data;
+ additional_data.mapping_update_flags = update_values;
+ additional_data.mapping_update_flags_inner_faces = update_values;
+ additional_data.mapping_update_flags_boundary_faces = update_values;
+ additional_data.mapping_update_flags_faces_by_cells = update_values;
+ additional_data.hold_all_faces_to_owned_cells = true;
+
+ MatrixFree<dim, Number, VectorizedArrayType> 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);
+
+ this->setup_vector(matrix_free);
+
+ matrix_free.cell_loop(&Test::cell_wise_operation, this, dst, src);
+
+ deallog << std::endl;
+ }
+
+private:
+ void
+ setup_vector(const MatrixFree<dim, Number, VectorizedArrayType> &data)
+ {
+ unsigned int n_cells = data.n_cell_batches() + data.n_ghost_cell_batches();
+ cell_ids.resize(n_cells);
+
+ for (unsigned int cell = 0; cell < n_cells; cell++)
+ for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell);
+ lane++)
+ cell_ids[cell][lane] = data.get_cell_iterator(cell, lane)->id();
+ }
+
+ void
+ cell_wise_operation(const MatrixFree<dim, Number, VectorizedArrayType> &data,
+ VectorType &,
+ const VectorType &,
+ const std::pair<unsigned int, unsigned int> &pair) const
+ {
+ FEEvaluation<dim, 1> fe_eval(data);
+ FEFaceEvaluation<dim, 1> fe_eval_m(data, true);
+ FEFaceEvaluation<dim, 1> fe_eval_p(data, false);
+
+ for (auto cell = pair.first; cell < pair.second; cell++)
+ {
+ for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell);
+ lane++)
+ {
+ fe_eval.reinit(cell);
+ const auto cell_data = fe_eval.read_cell_data(cell_ids);
+ Assert(cell_data[lane] == data.get_cell_iterator(cell, lane)->id(),
+ ExcInternalError());
+ // deallog << "c " << cell_data[lane] << std::endl;
+
+ for (unsigned int face = 0;
+ face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ {
+ fe_eval_m.reinit(cell, face);
+ fe_eval_p.reinit(cell, face);
+
+ const auto cell_data_m = fe_eval_m.read_cell_data(cell_ids);
+ Assert(cell_data_m[lane] ==
+ data.get_cell_iterator(cell, lane)->id(),
+ ExcInternalError());
+ // deallog << "m" << cell_data_m[lane] << std::endl;
+
+ const auto bids =
+ data.get_faces_by_cells_boundary_id(cell, face);
+
+ if (bids[lane] != numbers::invalid_unsigned_int)
+ continue;
+
+ const auto cell_data_p = fe_eval_p.read_cell_data(cell_ids);
+
+ Assert(
+ cell_data_p[lane] ==
+ data.get_cell_iterator(cell, lane)->neighbor(face)->id(),
+ ExcInternalError());
+ // deallog << "p " << cell_data_p[lane] << std::endl;
+ }
+ }
+ }
+ }
+
+ AlignedVector<std::array<CellId, VectorizedArrayType::n_array_elements>>
+ cell_ids;
+};
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+ MPILogInitAll all;
+
+ {
+ deallog.push("2d");
+ Test<2> runner;
+ runner.run(1);
+ deallog.pop();
+ }
+ {
+ deallog.push("3d");
+ Test<3> runner;
+ runner.run(1);
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:0:2d::
+DEAL:0:3d::
+
+DEAL:1:2d::
+DEAL:1:3d::
+
projects / dealii.git / commitdiff