From: Johannes Heinz Date: Mon, 27 Nov 2023 08:55:52 +0000 (+0100) Subject: Add FERemoteEvaluation X-Git-Tag: relicensing~217^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d79b9acb470d5fbfdb420b72f77c5631f527add8;p=dealii.git Add FERemoteEvaluation --- diff --git a/doc/news/changes/major/20231127HeinzMunch b/doc/news/changes/major/20231127HeinzMunch new file mode 100644 index 0000000000..2fadf547d2 --- /dev/null +++ b/doc/news/changes/major/20231127HeinzMunch @@ -0,0 +1,8 @@ +New: Introduce FERemoteEvaluation. This new class provides similar +interfaces as FEEvaluation, FEFaceEvaluation, and FEPointEvaluation. +It can be used to access values and/or gradients in quadrature points +at remote triangulations. In particular the class is designed such that +it can also be used for Nitsche-type mortaring or L2-Projection with +CGAL in the future. +
+(Johannes Heinz, Peter Munch, 2023/11/27) diff --git a/include/deal.II/matrix_free/fe_remote_evaluation.h b/include/deal.II/matrix_free/fe_remote_evaluation.h new file mode 100644 index 0000000000..94ef418901 --- /dev/null +++ b/include/deal.II/matrix_free/fe_remote_evaluation.h @@ -0,0 +1,1037 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_fe_remote_evaluation_h +#define dealii_matrix_free_fe_remote_evaluation_h + +#include + +#include +#include +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + /** + * A class that stores values and/or gradients, e.g. at quadrature points that + * are accessed inside matrix-free loops. + */ + template + struct PrecomputedEvaluationData + { + using value_type = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits::value_type; + + using gradient_type = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits::gradient_type; + + /** + * Values at quadrature points. + */ + AlignedVector values; + + /** + * Gradients at quadrature points. + */ + AlignedVector gradients; + }; + + /** + * A class that stores a CRS like structure to access + * PrecomputedEvaluationData. In most cases a one-level CRS structure + * is enough. In this case @c ptrs has to be constructed and the shift can be + * obtained with `get_shift(index)`. The field @c ptrs_ptrs stays empty. + * It is only filled if a two level structure is needed. In this case + * `get_shift(cell_index, face_number)` return the correct shift. + */ + struct PrecomputedEvaluationDataView + { + /** + * Get a pointer to data at @p index. + */ + unsigned int + get_shift(const unsigned int index) const; + + /** + * Get a pointer to data at (@p cell_index, @p face_number). + */ + unsigned int + get_shift(const unsigned int cell_index, + const unsigned int face_number) const; + + /** + * Get the number of stored values. + */ + unsigned int + size() const; + + /** + * This parameter can be used if indices do not start with 0. + */ + unsigned int start = 0; + + /** + * Pointers to @c ptrs. + */ + std::vector ptrs_ptrs; + + /** + * Pointers to data. + */ + std::vector ptrs; + }; + + /** + * A class helps to access PecomputedEvaluationData in a thread-safe + * manner. Each thread has to create this wrapper class on its own to + * avoid race-conditions during reinit(). + */ + template + class PrecomputedEvaluationDataAccessor + { + using value_type = + typename PrecomputedEvaluationData:: + value_type; + using gradient_type = + typename PrecomputedEvaluationData:: + gradient_type; + + public: + /** + * Constructor. + */ + PrecomputedEvaluationDataAccessor( + const PrecomputedEvaluationData &data, + const PrecomputedEvaluationDataView &view); + + /** + * Get the value at quadrature point @p q. The entity on which the values + * are defined is set via `reinit()`. + * + * @param[in] q Quadrature point at which the value is queried. + */ + const value_type + get_value(const unsigned int q) const; + + /** + * Get the gradients at quadrature point @p q. The entity on which the + * gradients are defined is set via `reinit()`. + * + * @param[in] q Quadrature point at which the gradient is queried. + */ + const gradient_type + get_gradient(const unsigned int q) const; + + /** + * This function has to be called before `get_value()` and/or + * `get_gradient()`. + * + * @param[in] index Entity index at which quadrature points are accessed. + * This can be, e.g., a cell index, a cell batch index, or a face batch + * index. + */ + void + reinit(const unsigned int index); + + /** + * Set cell and face number at which quadrature points are accessed. + * + * @param[in] index_0 cell index. + * @param[in] index_1 cell-local face number. + */ + void + reinit(const unsigned int index_0, const unsigned int index_1); + + private: + /** + * PrecomputedEvaluationDataView provides information where values are + * located. + */ + const PrecomputedEvaluationDataView &view; + + /** + * PrecomputedEvaluationData stores the actual values. + */ + const PrecomputedEvaluationData &data; + + /** + * Offset to data after last call of `reinit()`. + */ + unsigned int data_offset; + }; + +} // namespace internal + + + +/** + * Communication objects know about the communication pattern. + * In case of (matrix-free) cells batches or faces batches + * a RemotePointEvaluation object stores the location of the + * remote points. @c batch_id_n_entities relates these points + * to the corresponding quadrature points of entity batches. + * For this the field stores batch IDs and the number of entities + * in the batch. + */ +template +struct FERemoteCommunicationObjectEntityBatches +{ + /** + * Object that is reinitialized with the remote points we want to access. + */ + std::shared_ptr> rpe; + + /** + * A vector that stores the batch IDs and the number of active entries in the + * batch related to the the remote points. + */ + std::vector> batch_id_n_entities; + + /** + * Function that gives access to @c batch_id_n_entities. The function is only + * used for a unified access to the vector that stores additional information + * on the points that are processed via @c rpe. + */ + std::vector> + get_communication_object_pntrs() const; +}; + +/** + * Similar as @c FERemoteCommunicationObjectEntityBatches. + * To relate the points from @c RemotePointEvaluation to + * quadrature points on corresponding cells, cell iterators + * have to be stored. + */ +template +struct FERemoteCommunicationObjectCells +{ + /** + * Object that is reinitialized with the remote points we want to access. + */ + std::shared_ptr> rpe; + + /** + * A vector that stores the cell iterators related to the the remote points. + */ + std::vector::cell_iterator> cells; + + /** + * Function that gives access to @c cells. The function is only used + * for a unified access to the vector that stores additional information on + * the points that are processed via @c rpe. + */ + std::vector::cell_iterator> + get_communication_object_pntrs() const; +}; + +/** + * Similar as @c FERemoteCommunicationObjectCells. + * To relate the points from @c RemotePointEvaluation to + * quadrature points on corresponding faces, cell iterators + * and face numbers have to be stored. + */ +template +struct FERemoteCommunicationObjectFaces +{ + /** + * Object that is reinitialized with the remote points we want to access. + */ + std::shared_ptr> rpe; + + /** + * A vector that stores the cell iterators and the face number related to the + * the remote points. + */ + std::vector< + std::pair::cell_iterator, unsigned int>> + cell_face_nos; + + /** + * Function that gives access to @c cell_face_nos. The function is only used + * for a unified access to the vector that stores additional information on + * the points that are processed via @c rpe. + */ + std::vector< + std::pair::cell_iterator, unsigned int>> + get_communication_object_pntrs() const; +}; + +/** + * A class to fill the fields in PrecomputedEvaluationData. + */ +template +class FERemoteEvaluationCommunicator : public Subscriptor +{ +public: + /** + * This function stores given communication objects + * and constructs a @c PrecomputedEvaluationDataView + * object if remote points are related to matrix-free face batches. + */ + void + reinit_faces(const std::vector> + &comm_objects, + const std::pair &face_batch_range, + const std::vector> &quadrature_vector); + + /** + * This function stores given communication objects + * and constructs a @c PrecomputedEvaluationDataView + * object if remote points are related to faces of given + * cells. + */ + template + void + reinit_faces( + const std::vector> &comm_objects, + const IteratorRange &cell_iterator_range, + const std::vector>> &quadrature_vector); + + /** + * Fill the fields stored in PrecomputedEvaluationData. + */ + template + void + update_ghost_values( + PrecomputedEvaluationDataType &dst, + const MeshType &mesh, + const VectorType &src, + const EvaluationFlags::EvaluationFlags eval_flags, + const unsigned int first_selected_component, + const VectorTools::EvaluationFlags::EvaluationFlags vec_flags) const; + + /** + * Provide access to @c PrecomputedEvaluationDataView. + */ + const internal::PrecomputedEvaluationDataView & + get_view() const; + +private: + /** + * CRS like data structure that describes the data positions at given + * indices. + */ + internal::PrecomputedEvaluationDataView view; + + /** + * A variant for all possible communication objects. + */ + std::vector, + FERemoteCommunicationObjectCells, + FERemoteCommunicationObjectFaces>> + communication_objects; + + /** + * Functions in this class only deals with copying data to @c + * PrecomputedEvaluationData. + */ + class CopyInstructions + { + public: + /** + * Copy data from @c src to @c dst. Overload for @c + * FERemoteCommunicationObjectCells. + */ + template + static void + copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector &dst, + const std::vector &src, + const std::vector::cell_iterator> &cells); + + /** + * Copy data from @c src to @c dst. Overload for @c + * FERemoteCommunicationObjectFaces. + */ + template + static void + copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector &dst, + const std::vector &src, + const std::vector::cell_iterator, + unsigned int>> &cell_face_nos); + + /** + * Copy data from @c src to @c dst. Overload for @c + * FERemoteCommunicationObjectEntityBatches. + */ + template + static void + copy_data(const internal::PrecomputedEvaluationDataView &view, + AlignedVector &dst, + const std::vector &src, + const std::vector> + &batch_id_n_entities); + + private: + /** + * Copy data to the correct position in a @c VectorizedArray. + */ + template + static void + copy_data_entries(VectorizedArray &dst, + const unsigned int v, + const T1 &src); + + /** + * Similar as @c copy_data_entries() above. + */ + template + static void + copy_data_entries(Tensor> &dst, + const unsigned int v, + const Tensor &src); + + /** + * Similar as @c copy_data_entries() above. + */ + template + static void + copy_data_entries( + Tensor>> &dst, + const unsigned int v, + const Tensor> &src); + + /** + * Throw a runtime exception if @c copy_data_entries() has not been + * implemented for a given type. + */ + template + static void + copy_data_entries(T1 &, const unsigned int, const T2 &); + }; +}; + + + +/** + * Class to access data in a matrix-free loop for non-matching discretizations. + * Interfaces are named with FEEvaluation in mind. + * The main difference is, that `gather_evaluate()` updates and caches + * all values at once. Therefore, it has to be called on one thread before a + * matrix-free loop. + * + * To access values and gradients in a thread safe way, @c get_data_accessor() + * has to be called on every thread. It provides the functions `get_value()` and + * `get_gradient()`. + */ +template +class FERemoteEvaluation +{ +public: + /** + * The constructor needs a corresponding FERemoteEvaluationCommunicator + * which has to be setup outside of this class. This design choice is + * motivated since the same FERemoteEvaluationCommunicator can be used + * for different MeshTypes and number of components. + * + * @param[in] comm FERemoteEvaluationCommunicator. + * @param[in] mesh Triangulation or DoFHandler. + * @param[in] evaluation_flags Specify treatment of values at points which are + * found on multiple cells. + * @param[in] first_selected_component Select first component of evaluation in + * DoFHandlers with multiple components. + */ + template + FERemoteEvaluation(const FERemoteEvaluationCommunicator &comm, + const MeshType &mesh, + const unsigned int first_selected_component = 0, + const VectorTools::EvaluationFlags::EvaluationFlags + evaluation_flags = VectorTools::EvaluationFlags::avg); + + /** + * Update the data which can be accessed via `get_value()` and + * `get_gradient()`. + * + * @param[in] src Solution vector used to update data. + * @param[in] flags Evaluation flags. Currently supported are + * EvaluationFlags::values and EvaluationFlags::gradients. + */ + template + void + gather_evaluate(const VectorType &src, + const EvaluationFlags::EvaluationFlags flags); + + /** + * @c FERemoteEvaluation does not provide the functions `get_value()` and + *`get_gradient()`. To access values and/or gradients call @c + *get_data_accessor() on every thread, e.g., `auto remote_evaluator = + *get_data_accessor();` The returned object can be used as follows. + *@code + * for (unsigned int entity = range.first; entity < range.second; ++entity) + * { + * remote_evaluator.reinit(entity); + * for(unsigned int q : quadrature_point_indices()) + * remote_evaluator.get_value(q) + * } + *@endcode + */ + internal::PrecomputedEvaluationDataAccessor + get_data_accessor() const; + +private: + /** + * Use Triangulation as MeshType. + */ + void + set_mesh(const Triangulation &tria); + + /** + * Use DoFHandler as MeshType. + */ + void + set_mesh(const DoFHandler &dof_handler); + + /** + * Precomputed values and/or gradients at remote locations. + */ + internal::PrecomputedEvaluationData data; + + /** + * Underlying communicator which handles update of the ghost values. + */ + SmartPointer> comm; + + /** + * Pointer to MeshType if used with Triangulation. + */ + SmartPointer> tria; + + /** + * Pointer to MeshType if used with DoFHandler. + */ + SmartPointer> dof_handler; + + /** + * First selected component. + */ + const unsigned int first_selected_component; + + /** + * Flags that indicate which ghost values are updated. + */ + const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags; +}; + + + +namespace internal +{ + unsigned int + PrecomputedEvaluationDataView::get_shift(const unsigned int index) const + { + Assert(ptrs_ptrs.size() == 0, ExcMessage("Two level CRS set up")); + + Assert(index != numbers::invalid_unsigned_int, + ExcMessage("Index has to be valid!")); + + Assert(start <= index, ExcInternalError()); + AssertIndexRange(index - start, ptrs.size()); + return ptrs[index - start]; + } + + unsigned int + PrecomputedEvaluationDataView::get_shift(const unsigned int cell_index, + const unsigned int face_number) const + { + Assert(ptrs_ptrs.size() > 0, ExcMessage("No two level CRS set up")); + + Assert(cell_index != numbers::invalid_unsigned_int, + ExcMessage("Cell index has to be valid!")); + Assert(face_number != numbers::invalid_unsigned_int, + ExcMessage("Face number has to be valid!")); + + Assert(start <= cell_index, ExcInternalError()); + + AssertIndexRange(cell_index - start, ptrs_ptrs.size()); + const unsigned int face_index = ptrs_ptrs[cell_index - start] + face_number; + AssertIndexRange(face_index, ptrs.size()); + return ptrs[face_index]; + } + + unsigned int + PrecomputedEvaluationDataView::size() const + { + Assert(ptrs.size() > 0, ExcInternalError()); + return ptrs.back(); + } + + template + PrecomputedEvaluationDataAccessor:: + PrecomputedEvaluationDataAccessor( + const PrecomputedEvaluationData &data, + const PrecomputedEvaluationDataView &view) + : view(view) + , data(data) + , data_offset(numbers::invalid_unsigned_int) + {} + + template + const typename PrecomputedEvaluationData::value_type + PrecomputedEvaluationDataAccessor::get_value( + const unsigned int q) const + { + Assert(data_offset != numbers::invalid_unsigned_int, + ExcMessage("reinit() not called.")); + AssertIndexRange(data_offset + q, data.values.size()); + return data.values[data_offset + q]; + } + + template + const typename PrecomputedEvaluationData:: + gradient_type + PrecomputedEvaluationDataAccessor:: + get_gradient(const unsigned int q) const + { + Assert(data_offset != numbers::invalid_unsigned_int, + ExcMessage("reinit() not called.")); + AssertIndexRange(data_offset + q, data.gradients.size()); + return data.gradients[data_offset + q]; + } + + template + void + PrecomputedEvaluationDataAccessor::reinit( + const unsigned int index) + { + data_offset = view.get_shift(index); + } + + template + void + PrecomputedEvaluationDataAccessor::reinit( + const unsigned int index_0, + const unsigned int index_1) + { + data_offset = view.get_shift(index_0, index_1); + } + +} // namespace internal + + + +template +std::vector> +FERemoteCommunicationObjectEntityBatches::get_communication_object_pntrs() + const +{ + return batch_id_n_entities; +} + +template +std::vector::cell_iterator> +FERemoteCommunicationObjectCells::get_communication_object_pntrs() const +{ + return cells; +} + +template +std::vector::cell_iterator, unsigned int>> +FERemoteCommunicationObjectFaces::get_communication_object_pntrs() const +{ + return cell_face_nos; +} + + + +template +void +FERemoteEvaluationCommunicator::reinit_faces( + const std::vector> + &comm_objects, + const std::pair &face_batch_range, + const std::vector> &quadrature_vector) +{ + // erase type by converting to the base object + communication_objects.clear(); + for (const auto &co : comm_objects) + communication_objects.push_back(co); + + // fetch points and update communication patterns + const unsigned int n_cells = quadrature_vector.size(); + AssertDimension(n_cells, face_batch_range.second - face_batch_range.first); + + // construct view: + view.start = face_batch_range.first; + + view.ptrs.resize(n_cells + 1); + + view.ptrs[0] = 0; + for (unsigned int face = 0; face < n_cells; ++face) + { + view.ptrs[face + 1] = view.ptrs[face] + quadrature_vector[face].size(); + } +} + +template +template +void +FERemoteEvaluationCommunicator::reinit_faces( + const std::vector> &comm_objects, + const IteratorRange &cell_iterator_range, + const std::vector>> &quadrature_vector) +{ + // erase type + communication_objects.clear(); + for (const auto &co : comm_objects) + communication_objects.push_back(co); + + const unsigned int n_cells = quadrature_vector.size(); + AssertDimension(n_cells, + std::distance(cell_iterator_range.begin(), + cell_iterator_range.end())); + + // construct view: + auto &cell_ptrs = view.ptrs_ptrs; + auto &face_ptrs = view.ptrs; + + view.start = 0; + cell_ptrs.resize(n_cells); + unsigned int n_faces = 0; + for (const auto &cell : cell_iterator_range) + { + cell_ptrs[cell->active_cell_index()] = n_faces; + n_faces += cell->n_faces(); + } + + face_ptrs.resize(n_faces + 1); + face_ptrs[0] = 0; + for (const auto &cell : cell_iterator_range) + { + for (const auto &f : cell->face_indices()) + { + const unsigned int face_index = + cell_ptrs[cell->active_cell_index()] + f; + + face_ptrs[face_index + 1] = + face_ptrs[face_index] + + quadrature_vector[cell->active_cell_index()][f].size(); + } + } +} + +template +template +void +FERemoteEvaluationCommunicator::update_ghost_values( + PrecomputedEvaluationDataType &dst, + const MeshType &mesh, + const VectorType &src, + const EvaluationFlags::EvaluationFlags eval_flags, + const unsigned int first_selected_component, + const VectorTools::EvaluationFlags::EvaluationFlags vec_flags) const +{ + const bool has_ghost_elements = src.has_ghost_elements(); + + if (has_ghost_elements == false) + src.update_ghost_values(); + + + for (const auto &communication_object : communication_objects) + { + if (eval_flags & EvaluationFlags::values) + { + std::visit( + [&](const auto &obj) { + CopyInstructions::copy_data( + view, + dst.values, + VectorTools::point_values( + *obj.rpe, mesh, src, vec_flags, first_selected_component), + obj.get_communication_object_pntrs()); + }, + communication_object); + } + + if (eval_flags & EvaluationFlags::gradients) + { + std::visit( + [&](const auto &obj) { + CopyInstructions::copy_data( + view, + dst.gradients, + VectorTools::point_gradients( + *obj.rpe, mesh, src, vec_flags, first_selected_component), + obj.get_communication_object_pntrs()); + }, + communication_object); + } + + Assert(!(eval_flags & EvaluationFlags::hessians), ExcNotImplemented()); + } + + if (has_ghost_elements == false) + src.zero_out_ghost_values(); +} + +template +const internal::PrecomputedEvaluationDataView & +FERemoteEvaluationCommunicator::get_view() const +{ + return view; +} + +template +template +void +FERemoteEvaluationCommunicator::CopyInstructions::copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector &dst, + const std::vector &src, + const std::vector::cell_iterator> &cells) +{ + dst.resize(view.size()); + + unsigned int c = 0; + for (const auto &cell : cells) + { + for (unsigned int j = view.get_shift(cell->active_cell_index()); + j < view.get_shift(cell->active_cell_index() + 1); + ++j, ++c) + { + AssertIndexRange(j, dst.size()); + AssertIndexRange(c, src.size()); + dst[j] = src[c]; + } + } +} + +template +template +void +FERemoteEvaluationCommunicator::CopyInstructions::copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector &dst, + const std::vector &src, + const std::vector::cell_iterator, + unsigned int>> &cell_face_nos) +{ + dst.resize(view.size()); + + unsigned int c = 0; + for (const auto &[cell, f] : cell_face_nos) + { + for (unsigned int j = view.get_shift(cell->active_cell_index(), f); + j < view.get_shift(cell->active_cell_index(), f + 1); + ++j, ++c) + { + AssertIndexRange(j, dst.size()); + AssertIndexRange(c, src.size()); + + dst[j] = src[c]; + } + } +} + +template +template +void +FERemoteEvaluationCommunicator::CopyInstructions::copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector &dst, + const std::vector &src, + const std::vector> &batch_id_n_entities) +{ + dst.resize(view.size()); + + unsigned int c = 0; + for (const auto &[batch_id, n_entries] : batch_id_n_entities) + { + for (unsigned int v = 0; v < n_entries; ++v) + for (unsigned int j = view.get_shift(batch_id); + j < view.get_shift(batch_id + 1); + ++j, ++c) + { + AssertIndexRange(j, dst.size()); + AssertIndexRange(c, src.size()); + + copy_data_entries(dst[j], v, src[c]); + } + } +} + +template +template +void +FERemoteEvaluationCommunicator::CopyInstructions::copy_data_entries( + VectorizedArray &dst, + const unsigned int v, + const T1 &src) +{ + AssertIndexRange(v, n_lanes); + + dst[v] = src; +} + +template +template +void +FERemoteEvaluationCommunicator::CopyInstructions::copy_data_entries( + Tensor> &dst, + const unsigned int v, + const Tensor &src) +{ + AssertIndexRange(v, n_lanes); + + if constexpr (rank_ == 1) + { + for (unsigned int i = 0; i < dim_; ++i) + dst[i][v] = src[i]; + } + else + { + for (unsigned int i = 0; i < rank_; ++i) + for (unsigned int j = 0; j < dim_; ++j) + dst[i][j][v] = src[i][j]; + } +} + +template +template +void +FERemoteEvaluationCommunicator::CopyInstructions::copy_data_entries( + Tensor>> &dst, + const unsigned int v, + const Tensor> &src) +{ + if constexpr (rank_ == 1) + { + for (unsigned int i = 0; i < n_components_; ++i) + copy_data(dst[i], v, src[i]); + } + else + { + for (unsigned int i = 0; i < rank_; ++i) + for (unsigned int j = 0; j < n_components_; ++j) + dst[i][j][v] = src[i][j]; + } +} + +template +template +void +FERemoteEvaluationCommunicator::CopyInstructions::copy_data_entries( + T1 &, + const unsigned int, + const T2 &) +{ + Assert(false, + ExcMessage( + "copy_data_entries() not implemented for given arguments.")); +} + + + +template +template +FERemoteEvaluation::FERemoteEvaluation( + const FERemoteEvaluationCommunicator &comm, + const MeshType &mesh, + const unsigned int first_selected_component, + const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags) + : comm(&comm) + , first_selected_component(first_selected_component) + , evaluation_flags(evaluation_flags) +{ + set_mesh(mesh); +} + +template +template +void +FERemoteEvaluation::gather_evaluate( + const VectorType &src, + const EvaluationFlags::EvaluationFlags flags) +{ + if (tria) + { + Assert(n_components == 1, ExcNotImplemented()); + comm->template update_ghost_values(this->data, + *tria, + src, + flags, + first_selected_component, + evaluation_flags); + } + else if (dof_handler) + { + comm->template update_ghost_values(this->data, + *dof_handler, + src, + flags, + first_selected_component, + evaluation_flags); + } + else + Assert(false, ExcNotImplemented()); +} + +template +internal::PrecomputedEvaluationDataAccessor +FERemoteEvaluation::get_data_accessor() const +{ + internal::PrecomputedEvaluationDataAccessor data_accessor(data, + comm->get_view()); + return data_accessor; +} + +template +void +FERemoteEvaluation::set_mesh( + const Triangulation &tria) +{ + this->tria = &tria; +} + +template +void +FERemoteEvaluation::set_mesh( + const DoFHandler &dof_handler) +{ + this->dof_handler = &dof_handler; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/matrix_free/fe_remote_evaluation_01.cc b/tests/matrix_free/fe_remote_evaluation_01.cc new file mode 100644 index 0000000000..278d1a6a02 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_01.cc @@ -0,0 +1,424 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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. +// +// --------------------------------------------------------------------- + +// This test checks the access of values on boundary faces using +// FERemoteEvaluation. + +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +using Number = double; + +using VectorType = LinearAlgebra::distributed::Vector; + +constexpr unsigned int dim = 2; + +template +struct Coordinates : public Function +{ + Coordinates(unsigned int offset) + : Function(dim + offset, 0.0) + , offset(offset) + {} + + // Function that returns the coordinates of given point. + // In case an offset is provided, 0.0 is returned for values + // below the offset and the coordinates are returned for + // components >= offset. This is used to check the functionality + // for first selected component. + double + value(const Point &p, const unsigned int component) const final + { + if (component < offset) + return 0.0; + return p[component - offset]; + } + + unsigned int offset; +}; + +FERemoteEvaluationCommunicator +construct_comm_for_face_batches(const MatrixFree &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectEntityBatches co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // face + std::vector> quadrature_vector( + matrix_free.n_boundary_face_batches()); + + // Points that are searched by rpe. + std::vector> points; + + FEFaceEvaluation phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + + phi.reinit(face); + + const unsigned int n_faces = + matrix_free.n_active_entries_per_face_batch(face); + + co.batch_id_n_entities.push_back(std::make_pair(face, n_faces)); + + // Append the quadrature points to the points we need to search + // for. + for (unsigned int v = 0; v < n_faces; ++v) + { + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + + points.push_back(temp); + } + } + + // append quadrature of correct size + quadrature_vector[bface] = Quadrature(phi.n_q_points); + } + + // use rpe to search for stored points + auto rpe = + std::make_shared>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator remote_communicator; + remote_communicator.reinit_faces({co}, face_batch_range, quadrature_vector); + + return remote_communicator; +} + +void +test_face_batches(const MatrixFree &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_face_batches(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + Coordinates(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation> phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::values); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + phi_remote.reinit(face); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto error = + phi.quadrature_point(q) - phi_remote.get_value(q); + + // only consider active entries in global error + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + AssertThrow(std::abs(error[d][v]) < 1e-13, + ExcMessage("Error too large.")); + } + } + }; + + matrix_free.template loop( + {}, {}, boundary_function, dst, src, true); + } +} + + +FERemoteEvaluationCommunicator +construct_comm_for_cell_face_nos(const MatrixFree &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectFaces co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // cell, face pair + std::vector>> quadrature_vector; + for (const auto &cell : matrix_free.get_dof_handler() + .get_triangulation() + .active_cell_iterators()) + quadrature_vector.emplace_back( + std::vector>(cell->n_faces())); + + // Points that are searched by rpe. + std::vector> points; + + FEFaceEvaluation phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + phi.reinit(face); + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = matrix_free.get_face_iterator(face, v, true); + co.cell_face_nos.push_back(std::make_pair(cell, f)); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + points.push_back(temp); + } + + // append quadrature of correct size + quadrature_vector[cell->active_cell_index()][f] = + Quadrature(phi.n_q_points); + } + } + + // use rpe to search for stored points + auto rpe = + std::make_shared>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator remote_communicator; + remote_communicator.reinit_faces( + {co}, + matrix_free.get_dof_handler().get_triangulation().active_cell_iterators(), + quadrature_vector); + + return remote_communicator; +} + +void +test_cell_face_nos(const MatrixFree &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_cell_face_nos(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + Coordinates(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::values); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = + matrix_free.get_face_iterator(face, v, true); + phi_remote.reinit(cell->active_cell_index(), f); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Tensor<1, dim, Number> temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + + const auto error = temp - phi_remote.get_value(q); + for (unsigned int d = 0; d < dim; ++d) + AssertThrow(std::abs(error[d]) < 1e-13, + ExcMessage("Error too large.")); + } + } + } + }; + + matrix_free.template loop( + {}, {}, boundary_function, dst, src, true); + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + MPILogInitAll all; + + // 1) setup MatrixFree for different finite elements, i.e. + // FESystem(FE_DGQ, dim + first_selected_component). + + const unsigned int refinements = 2; + const unsigned int degree = 5; + + // Construct triangulation. + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(refinements); + + // Create DoFHandlers with different number of components to + // check different first_selected components + std::vector>> dof_handlers; + unsigned int max_first_selected_comp = 3; + for (unsigned int i = 0; i < max_first_selected_comp; ++i) + { + dof_handlers.push_back(std::make_unique>(tria)); + dof_handlers.back()->distribute_dofs( + FESystem(FE_DGQ(degree), dim + i)); + } + + // Constraints are not considered in this test + AffineConstraints constraints; + constraints.close(); + + // Setup MatrixFree + typename MatrixFree::AdditionalData data; + data.mapping_update_flags = update_values; + // We are not considering inner faces in this test. Therefore, + // data.mapping_update_flags_inner_faces is not set. + data.mapping_update_flags_boundary_faces = + update_quadrature_points | update_values; + + + // Construct vector of raw pointers needed to init matrix free + std::vector *> dof_handlers_mf; + std::vector *> constraints_mf; + for (auto &dh : dof_handlers) + { + dof_handlers_mf.push_back(dh.get()); + constraints_mf.push_back(&constraints); + } + + MatrixFree matrix_free; + matrix_free.reinit(MappingQ1(), + dof_handlers_mf, + constraints_mf, + QGauss(degree + 1), + data); + + // 2) test to access remote values at points in face batches + test_face_batches(matrix_free, max_first_selected_comp); + + // 3) test to access remote values at points on faces of cells + test_cell_face_nos(matrix_free, max_first_selected_comp); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=1.output b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=1.output @@ -0,0 +1,2 @@ + +DEAL:0::OK diff --git a/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=2.output b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..52103938b0 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=2.output @@ -0,0 +1,5 @@ + +DEAL:0::OK + +DEAL:1::OK + diff --git a/tests/matrix_free/fe_remote_evaluation_02.cc b/tests/matrix_free/fe_remote_evaluation_02.cc new file mode 100644 index 0000000000..d4b5fdb14a --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_02.cc @@ -0,0 +1,428 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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. +// +// --------------------------------------------------------------------- + +// This test checks the access of gradients on boundary faces using +// FERemoteEvaluation. + +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +using Number = double; + +using VectorType = LinearAlgebra::distributed::Vector; + +constexpr unsigned int dim = 2; + +template +struct QuadraticFunction : public Function +{ + QuadraticFunction(unsigned int offset) + : Function(dim + offset, 0.0) + , offset(offset) + {} + + // Function that returns the coordinates of given point. + // In case an offset is provided, 0.0 is returned for values + // below the offset and the coordinates are returned for + // components >= offset. This is used to check the functionality + // for first selected component. + double + value(const Point &p, const unsigned int component) const final + { + if (component < offset) + return 0.0; + return p[component - offset] * p[component - offset]; + } + + unsigned int offset; +}; + +FERemoteEvaluationCommunicator +construct_comm_for_face_batches(const MatrixFree &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectEntityBatches co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // face + std::vector> quadrature_vector( + matrix_free.n_boundary_face_batches()); + + // Points that are searched by rpe. + std::vector> points; + + FEFaceEvaluation phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + + phi.reinit(face); + + const unsigned int n_faces = + matrix_free.n_active_entries_per_face_batch(face); + + co.batch_id_n_entities.push_back(std::make_pair(face, n_faces)); + + // Append the quadrature points to the points we need to search + // for. + for (unsigned int v = 0; v < n_faces; ++v) + { + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + + points.push_back(temp); + } + } + + // append quadrature of correct size + quadrature_vector[bface] = Quadrature(phi.n_q_points); + } + + // use rpe to search for stored points + auto rpe = + std::make_shared>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator remote_communicator; + remote_communicator.reinit_faces({co}, face_batch_range, quadrature_vector); + + return remote_communicator; +} + +void +test_face_batches(const MatrixFree &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_face_batches(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + QuadraticFunction(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation> phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::gradients); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + phi.gather_evaluate(src, EvaluationFlags::gradients); + phi_remote.reinit(face); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto error = + phi.get_gradient(q) - phi_remote.get_gradient(q); + + // only consider active entries in global error + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + AssertThrow(std::abs(error[i][d][v]) < 1e-13, + ExcMessage("Error too large.")); + } + } + }; + + matrix_free.template loop( + {}, {}, boundary_function, dst, src, true); + } +} + + +FERemoteEvaluationCommunicator +construct_comm_for_cell_face_nos(const MatrixFree &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectFaces co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // cell, face pair + std::vector>> quadrature_vector; + for (const auto &cell : matrix_free.get_dof_handler() + .get_triangulation() + .active_cell_iterators()) + quadrature_vector.emplace_back( + std::vector>(cell->n_faces())); + + // Points that are searched by rpe. + std::vector> points; + + FEFaceEvaluation phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + phi.reinit(face); + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = matrix_free.get_face_iterator(face, v, true); + co.cell_face_nos.push_back(std::make_pair(cell, f)); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + points.push_back(temp); + } + + // append quadrature of correct size + quadrature_vector[cell->active_cell_index()][f] = + Quadrature(phi.n_q_points); + } + } + + // use rpe to search for stored points + auto rpe = + std::make_shared>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator remote_communicator; + remote_communicator.reinit_faces( + {co}, + matrix_free.get_dof_handler().get_triangulation().active_cell_iterators(), + quadrature_vector); + + return remote_communicator; +} + +void +test_cell_face_nos(const MatrixFree &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_cell_face_nos(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + QuadraticFunction(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::gradients); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + phi.gather_evaluate(src, EvaluationFlags::gradients); + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = + matrix_free.get_face_iterator(face, v, true); + phi_remote.reinit(cell->active_cell_index(), f); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto grad = phi.get_gradient(q); + Tensor<2, dim, Number> temp; + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < dim; ++j) + temp[i][j] = grad[i][j][v]; + + const auto error = temp - phi_remote.get_gradient(q); + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < dim; ++j) + AssertThrow(std::abs(error[i][j]) < 1e-13, + ExcMessage("Error too large.")); + } + } + } + }; + + matrix_free.template loop( + {}, {}, boundary_function, dst, src, true); + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + MPILogInitAll all; + + // 1) setup MatrixFree for different finite elements, i.e. + // FESystem(FE_DGQ, dim + first_selected_component). + + const unsigned int refinements = 2; + const unsigned int degree = 5; + + // Construct triangulation. + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(refinements); + + // Create DoFHandlers with different number of components to + // check different first_selected components + std::vector>> dof_handlers; + unsigned int max_first_selected_comp = 3; + for (unsigned int i = 0; i < max_first_selected_comp; ++i) + { + dof_handlers.push_back(std::make_unique>(tria)); + dof_handlers.back()->distribute_dofs( + FESystem(FE_DGQ(degree), dim + i)); + } + + // Constraints are not considered in this test + AffineConstraints constraints; + constraints.close(); + + // Setup MatrixFree + typename MatrixFree::AdditionalData data; + data.mapping_update_flags = update_gradients; + // We are not considering inner faces in this test. Therefore, + // data.mapping_update_flags_inner_faces is not set. + data.mapping_update_flags_boundary_faces = + update_quadrature_points | update_gradients; + + + // Construct vector of raw pointers needed to init matrix free + std::vector *> dof_handlers_mf; + std::vector *> constraints_mf; + for (auto &dh : dof_handlers) + { + dof_handlers_mf.push_back(dh.get()); + constraints_mf.push_back(&constraints); + } + + MatrixFree matrix_free; + matrix_free.reinit(MappingQ1(), + dof_handlers_mf, + constraints_mf, + QGauss(degree + 1), + data); + + // 2) test to access remote values at points in face batches + test_face_batches(matrix_free, max_first_selected_comp); + + // 3) test to access remote values at points on faces of cells + test_cell_face_nos(matrix_free, max_first_selected_comp); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=1.output b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=1.output @@ -0,0 +1,2 @@ + +DEAL:0::OK diff --git a/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=2.output b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..52103938b0 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=2.output @@ -0,0 +1,5 @@ + +DEAL:0::OK + +DEAL:1::OK +