#include <deal.II/base/mpi_remote_point_evaluation.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/fe_point_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/numerics/vector_tools.h>
+#include <algorithm>
#include <variant>
DEAL_II_NAMESPACE_OPEN
+/**
+ * The namespace collects convenience functions to setup
+ * FERemoteEvaluationCommunicator for typical use cases.
+ */
+namespace Utilities
+{
+ /**
+ * A factory function for the FERemoteEvaluationCommunicator in the case of
+ * point-to-point interpolation.
+ *
+ * @param[in] matrix_free MatrixFree object that is used to distribute
+ * quadrature points at non-matching faces.
+ * @param[in] non_matching_faces_marked_vertices A vector of boundary face IDs
+ * that relate to non-matching interfaces. Each boundary face ID is
+ * accompanied by a lambda that marks the vertices of cells which are
+ * considered during the search of remote points (quadrature points of faces
+ * with given boundary face ID).
+ * @param[in] quad_no Quadrature number in @p matrix_free.
+ * @param[in] dof_no DoFHandler number in @p matrix_free.
+ * @param[in] tolerance Tolerance to find remote points.
+ */
+ template <int dim,
+ typename Number,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+ FERemoteEvaluationCommunicator<dim>
+ setup_remote_communicator_faces_point_to_point_interpolation(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const std::vector<
+ std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+ &non_matching_faces_marked_vertices,
+ const unsigned int quad_no = 0,
+ const unsigned int dof_no = 0,
+ const double tolerance = 1e-9);
+
+
+
+ /**
+ * A factory function for the FERemoteEvaluationCommunicator in the case of
+ * Nitsche-type mortaring.
+ *
+ * @param[in] matrix_free MatrixFree object that is used to distribute
+ * @param[in] non_matching_faces_marked_vertices A vector of boundary face IDs
+ * that relate to non-matching interfaces. Each boundary face ID is
+ * accompanied by a lambda that marks the vertices of cells which are
+ * considered during the process of computing intersections.
+ * @param[in] n_q_pnts_1D The number of 1D quadrature points per intersection.
+ * @param[in] dof_no DoFHandler number in @p matrix_free.
+ * distributed to each intersection (given for one coordinate direction).
+ * @param[in] nm_mapping_info In case nm_mapping_info is not a `nullptr` it is
+ * setup for cell based loops, such that it can be used in combination with
+ * FERemoteEvaluation.
+ * @param[in] tolerance Tolerance to find intersections.
+ */
+ template <int dim,
+ typename Number,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+ FERemoteEvaluationCommunicator<dim>
+ setup_remote_communicator_faces_nitsche_type_mortaring(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const std::vector<
+ std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+ &non_matching_faces_marked_vertices,
+ const unsigned int n_q_pnts_1D,
+ const unsigned int dof_no = 0,
+ NonMatching::MappingInfo<dim, dim, Number> *nm_mapping_info = nullptr,
+ const double tolerance = 1e-9);
+} // namespace Utilities
+
+
+
/**
* Class to access data in a matrix-free loop for non-matching discretizations.
* Interfaces are named with FEEvaluation in mind.
+namespace Utilities
+{
+ template <int dim, typename Number, typename VectorizedArrayType>
+ FERemoteEvaluationCommunicator<dim>
+ setup_remote_communicator_faces_point_to_point_interpolation(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const std::vector<
+ std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+ &non_matching_faces_marked_vertices,
+ const unsigned int quad_no,
+ const unsigned int dof_no,
+ const double tolerance)
+ {
+ const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
+ const auto &tria = dof_handler.get_triangulation();
+ const auto &mapping = *matrix_free.get_mapping_info().mapping;
+
+ // Communication objects know about the communication pattern. I.e.,
+ // they know about the cells and quadrature points that have to be
+ // evaluated at remote faces. This information is given via
+ // RemotePointEvaluation. Additionally, the communication objects
+ // have to be able to match the quadrature points of the remote
+ // points (that provide exterior information) to the quadrature points
+ // defined at the interior cell. In case of point-to-point interpolation
+ // a vector of pairs with face batch Ids and the number of faces in the
+ // batch is needed. @c FERemoteCommunicationObjectEntityBatches
+ // is a container to store this information.
+ //
+ // We need multiple communication objects (one for each non-matching face
+ // ID).
+ std::vector<FERemoteCommunicationObjectEntityBatches<dim>> comm_objects;
+
+ // Additionally to the communication objects we need a vector
+ // that stores quadrature rules for every face batch.
+ // The quadrature can be empty in case of non non-matching faces,
+ // i.e. boundary faces. Internally this information is needed to correctly
+ // access values over multiple communication objects.
+ std::vector<Quadrature<dim>> global_quadrature_vector(
+ matrix_free.n_boundary_face_batches());
+
+ // Get the range of face batches we have to look at during construction of
+ // the communication objects. We only have to look at boundary faces.
+ 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());
+
+ // Iterate over all non-matching face IDs.
+ for (const auto &[nm_face, marked_vertices] :
+ non_matching_faces_marked_vertices)
+ {
+ // Construct the communication object for every face ID:
+ // 1) RemotePointEvaluation with user specified function for marked
+ // vertices.
+ auto rpe = std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(
+ tolerance, false, 0, marked_vertices);
+
+ // 2) Face batch IDs and number of faces in batch.
+ std::vector<std::pair<unsigned int, unsigned int>>
+ face_batch_id_n_faces;
+
+ // Points that are searched by rpe.
+ std::vector<Point<dim>> points;
+
+ // Temporarily setup FEFaceEvaluation to access the quadrature points at
+ // the faces on the non-matching interface.
+ FEFaceEvaluation<dim, -1, 0, 1, Number> phi(matrix_free,
+ true,
+ dof_no,
+ quad_no);
+
+ // Iterate over the boundary faces.
+ for (unsigned int bface = 0;
+ bface < face_batch_range.second - face_batch_range.first;
+ ++bface)
+ {
+ const unsigned int face = face_batch_range.first + bface;
+
+ if (matrix_free.get_boundary_id(face) == nm_face)
+ {
+ phi.reinit(face);
+
+ // If @c face is on the current side of the non-matching
+ // interface. Add the face batch ID and the number of faces in
+ // the batch to the corresponding data structure.
+ const unsigned int n_faces =
+ matrix_free.n_active_entries_per_face_batch(face);
+ face_batch_id_n_faces.emplace_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<dim> temp;
+ for (unsigned int i = 0; i < dim; ++i)
+ temp[i] = point[i][v];
+
+ points.push_back(temp);
+ }
+ }
+
+ // Insert a quadrature rule of correct size into the global
+ // quadrature vector. First check that each face is only
+ // considered once.
+ Assert(global_quadrature_vector[bface].size() == 0,
+ ExcMessage(
+ "Quadrature for given face already provided."));
+
+ global_quadrature_vector[bface] =
+ Quadrature<dim>(phi.n_q_points);
+ }
+ }
+
+ // Reinit RPE and ensure all points are found.
+ rpe->reinit(points, tria, mapping);
+ Assert(rpe->all_points_found(),
+ ExcMessage("Not all remote points found."));
+
+ // Add communication object to the list of objects.
+ FERemoteCommunicationObjectEntityBatches<dim> co;
+ co.batch_id_n_entities = face_batch_id_n_faces;
+ co.rpe = rpe;
+ comm_objects.push_back(co);
+ }
+
+ // Reinit the communicator `FERemoteEvaluationCommunicator`
+ // with the communication objects.
+ FERemoteEvaluationCommunicator<dim> remote_communicator;
+
+ remote_communicator.reinit_faces(comm_objects,
+ face_batch_range,
+ global_quadrature_vector);
+
+ return remote_communicator;
+ }
+
+
+
+ template <int dim, typename Number, typename VectorizedArrayType>
+ FERemoteEvaluationCommunicator<dim>
+ setup_remote_communicator_faces_nitsche_type_mortaring(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const std::vector<
+ std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+ &non_matching_faces_marked_vertices,
+ const unsigned int n_q_pnts_1D,
+ const unsigned int dof_no,
+ NonMatching::MappingInfo<dim, dim, Number> *nm_mapping_info,
+ const double tolerance)
+ {
+ const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
+ const auto &tria = dof_handler.get_triangulation();
+ const auto &mapping = *matrix_free.get_mapping_info().mapping;
+
+ constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
+
+ std::pair<unsigned int, unsigned int> face_range =
+ std::make_pair(matrix_free.n_inner_face_batches(),
+ matrix_free.n_inner_face_batches() +
+ matrix_free.n_boundary_face_batches());
+
+ std::vector<Quadrature<dim - 1>> global_quadrature_vector(
+ (matrix_free.n_inner_face_batches() +
+ matrix_free.n_boundary_face_batches()) *
+ n_lanes);
+
+ // In case of Nitsche-type mortaring a vector of face indices is
+ // needed as communication object.
+ // @c FERemoteCommunicationObjectFaces is a container to store this
+ // information.
+ //
+ // We need multiple communication objects (one for each non-matching face
+ // ID).
+ std::vector<FERemoteCommunicationObject<dim>> comm_objects;
+
+ // Create bounding boxes and GridTools::Cache which is needed in
+ // the following loop.
+ std::vector<BoundingBox<dim>> local_boxes;
+ for (const auto &cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned())
+ local_boxes.emplace_back(mapping.get_bounding_box(cell));
+
+ // Create r-tree of bounding boxes
+ const auto local_tree = pack_rtree(local_boxes);
+
+ // Compress r-tree to a minimal set of bounding boxes
+ std::vector<std::vector<BoundingBox<dim>>> global_bboxes(1);
+ global_bboxes[0] = extract_rtree_level(local_tree, 0);
+
+ const GridTools::Cache<dim, dim> cache(tria, mapping);
+
+ // Iterate over all sides of the non-matching interface.
+ for (const auto &[nm_face, marked_vertices] :
+ non_matching_faces_marked_vertices)
+ {
+ // 1) compute cell face pairs
+ std::vector<
+ std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>
+ cell_face_pairs;
+
+ std::vector<unsigned int> indices;
+
+ for (unsigned int face = face_range.first; face < face_range.second;
+ ++face)
+ {
+ if (matrix_free.get_boundary_id(face) == nm_face)
+ {
+ for (unsigned int v = 0;
+ v < matrix_free.n_active_entries_per_face_batch(face);
+ ++v)
+ {
+ const auto &[c, f] = matrix_free.get_face_iterator(face, v);
+
+ cell_face_pairs.emplace_back(std::make_pair(c, f));
+ indices.push_back(face * n_lanes + v);
+ }
+ }
+ }
+
+ // 2) Create RPE.
+ // In the Nitsche-type case we do not collect points for the setup
+ // of RemotePointEvaluation. Instead we compute intersections between
+ // the faces and setup RemotePointEvaluation with the computed
+ // intersections.
+
+ // Build intersection requests. Intersection requests
+ // correspond to vertices at faces.
+ std::vector<std::vector<Point<dim>>> intersection_requests;
+ for (const auto &[cell, f] : cell_face_pairs)
+ {
+ std::vector<Point<dim>> vertices(cell->face(f)->n_vertices());
+ std::copy_n(mapping.get_vertices(cell, f).begin(),
+ cell->face(f)->n_vertices(),
+ vertices.begin());
+ intersection_requests.emplace_back(vertices);
+ }
+
+ // Compute intersection data with user specified function for marked
+ // vertices.
+ auto intersection_data =
+ GridTools::internal::distributed_compute_intersection_locations<
+ dim - 1>(cache,
+ intersection_requests,
+ global_bboxes,
+ marked_vertices(),
+ tolerance);
+
+ // Convert to RPE.
+ std::vector<Quadrature<dim>> mapped_quadratures_recv_comp;
+
+ auto rpe =
+ std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>();
+ rpe->reinit(
+ intersection_data
+ .template convert_to_distributed_compute_point_locations_internal<
+ dim>(n_q_pnts_1D, tria, mapping, &mapped_quadratures_recv_comp),
+ tria,
+ mapping);
+
+ // 3) Fill global quadrature vector.
+ for (unsigned int i = 0; i < intersection_requests.size(); ++i)
+ {
+ const auto idx = indices[i];
+
+ // We do not use a structural binding here, since with
+ // C++17 caputuring structural bindings in lambdas leads
+ // to an ill formed program.
+ const auto &cell = std::get<0>(cell_face_pairs[i]);
+ const auto &f = std::get<1>(cell_face_pairs[i]);
+
+ std::vector<Point<dim - 1>> q_points;
+ std::vector<double> weights;
+ for (unsigned int ptr = intersection_data.recv_ptrs[i];
+ ptr < intersection_data.recv_ptrs[i + 1];
+ ++ptr)
+ {
+ const auto &quad = mapped_quadratures_recv_comp[ptr];
+
+ const auto &ps = quad.get_points();
+ std::transform(
+ ps.begin(),
+ ps.end(),
+ std::back_inserter(q_points),
+ [&](const Point<dim> &p) {
+ return mapping.project_real_point_to_unit_point_on_face(
+ cell, f, p);
+ });
+
+ const auto &ws = quad.get_weights();
+ weights.insert(weights.end(), ws.begin(), ws.end());
+ }
+ Quadrature<dim - 1> quad(q_points, weights);
+
+ Assert(global_quadrature_vector[idx].size() == 0,
+ ExcMessage("Quadrature for given face already provided."));
+
+ global_quadrature_vector[idx] = quad;
+ }
+
+ // Add communication object.
+ FERemoteCommunicationObject<dim> co;
+ co.indices = indices;
+ co.rpe = rpe;
+ comm_objects.push_back(co);
+ }
+
+ // Reinit the communicator with the communication objects.
+ FERemoteEvaluationCommunicator<dim> remote_communicator;
+
+ remote_communicator.reinit_faces(
+ comm_objects,
+ std::make_pair(0, global_quadrature_vector.size()),
+ global_quadrature_vector);
+
+ if (nm_mapping_info != nullptr)
+ {
+ std::vector<
+ std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>>
+ vector_face_accessors;
+ vector_face_accessors.reserve((matrix_free.n_inner_face_batches() +
+ matrix_free.n_boundary_face_batches()) *
+ n_lanes);
+
+ // fill container for inner face batches
+ unsigned int face_batch = 0;
+ for (; face_batch < matrix_free.n_inner_face_batches(); ++face_batch)
+ {
+ for (unsigned int v = 0; v < n_lanes; ++v)
+ {
+ if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
+ vector_face_accessors.push_back(
+ matrix_free.get_face_iterator(face_batch, v));
+ else
+ vector_face_accessors.push_back(
+ matrix_free.get_face_iterator(face_batch, 0));
+ }
+ }
+ // and boundary face batches
+ for (; face_batch < (matrix_free.n_inner_face_batches() +
+ matrix_free.n_boundary_face_batches());
+ ++face_batch)
+ {
+ for (unsigned int v = 0; v < n_lanes; ++v)
+ {
+ if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
+ vector_face_accessors.push_back(
+ matrix_free.get_face_iterator(face_batch, v));
+ else
+ vector_face_accessors.push_back(
+ matrix_free.get_face_iterator(face_batch, 0));
+ }
+ }
+
+ nm_mapping_info->reinit_faces(vector_face_accessors,
+ global_quadrature_vector);
+ }
+
+ return remote_communicator;
+ }
+} // namespace Utilities
+
+
+
template <int dim, int n_components, typename value_type>
template <typename MeshType>
FERemoteEvaluation<dim, n_components, value_type>::FERemoteEvaluation(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 correct setup of FERemoteEvaluationCommunicator
+// via the factory functions
+// setup_remote_communicator_faces_point_to_point_interpolation()
+// and setup_remote_communicator_faces_nitsche_type_mortaring().
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/la_parallel_block_vector.h>
+
+#include <deal.II/matrix_free/fe_remote_evaluation.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+using Number = double;
+
+using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+constexpr unsigned int dim = 2;
+
+template <int dim>
+struct Coordinates : public Function<dim>
+{
+ Coordinates()
+ : Function<dim>(dim, 0.0)
+ {}
+
+ double
+ value(const Point<dim> &p, const unsigned int component) const final
+ {
+ if (component == 0)
+ return p[component];
+ return p[component];
+ }
+};
+
+
+void
+test_point_to_point(const MatrixFree<dim, Number> &matrix_free)
+{
+ const Triangulation<dim> &tria =
+ matrix_free.get_dof_handler().get_triangulation();
+
+ std::vector<std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+ non_matching_faces_marked_vertices;
+
+ non_matching_faces_marked_vertices.emplace_back(std::make_pair(1, [&]() {
+ std::vector<bool> mask(tria.n_vertices(), true);
+
+ for (const auto &cell : tria.active_cell_iterators())
+ for (auto const &f : cell->face_indices())
+ if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1)
+ for (const auto v : cell->vertex_indices())
+ mask[cell->vertex_index(v)] = false;
+
+ return mask;
+ }));
+
+ auto remote_communicator =
+ Utilities::setup_remote_communicator_faces_point_to_point_interpolation(
+ matrix_free, non_matching_faces_marked_vertices);
+
+ // 1) Allocate vectors and interpolate grid coordinates in DoFs.
+ VectorType src;
+ VectorType dst;
+ matrix_free.initialize_dof_vector(src);
+ matrix_free.initialize_dof_vector(dst);
+
+ VectorTools::interpolate(*matrix_free.get_mapping_info().mapping,
+ matrix_free.get_dof_handler(),
+ Coordinates<dim>(),
+ src);
+
+ // 2) Construct remote evaluator...
+ FERemoteEvaluation<dim, dim, VectorizedArray<Number>> phi_remote_eval(
+ remote_communicator, matrix_free.get_dof_handler());
+
+ // ...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<dim, -1, 0, dim, Number> phi(matrix_free, true);
+ auto phi_remote = phi_remote_eval.get_data_accessor();
+
+ for (unsigned int face = face_range.first; face < face_range.second;
+ ++face)
+ {
+ if (matrix_free.get_boundary_id(face) == 1)
+ {
+ 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<VectorType, VectorType>(
+ {}, {}, boundary_function, dst, src, true);
+}
+
+void
+test_nitsche_type_mortaring(const MatrixFree<dim, Number> &matrix_free)
+{
+ const Triangulation<dim> &tria =
+ matrix_free.get_dof_handler().get_triangulation();
+
+ std::vector<std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+ non_matching_faces_marked_vertices;
+
+ non_matching_faces_marked_vertices.emplace_back(std::make_pair(1, [&]() {
+ std::vector<bool> mask(tria.n_vertices(), true);
+
+ for (const auto &cell : tria.active_cell_iterators())
+ for (auto const &f : cell->face_indices())
+ if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1)
+ for (const auto v : cell->vertex_indices())
+ mask[cell->vertex_index(v)] = false;
+
+ return mask;
+ }));
+
+ typename NonMatching::MappingInfo<dim, dim, Number>::AdditionalData
+ additional_data;
+ additional_data.use_global_weights = true;
+
+ NonMatching::MappingInfo<dim, dim, Number> nm_mapping_info(
+ *matrix_free.get_mapping_info().mapping,
+ update_quadrature_points,
+ additional_data);
+
+ auto remote_communicator =
+ Utilities::setup_remote_communicator_faces_nitsche_type_mortaring(
+ matrix_free, non_matching_faces_marked_vertices, 4, 0, &nm_mapping_info);
+
+ // 1) Allocate vectors and interpolate grid coordinates in DoFs.
+ VectorType src;
+ VectorType dst;
+ matrix_free.initialize_dof_vector(src);
+ matrix_free.initialize_dof_vector(dst);
+
+ VectorTools::interpolate(*matrix_free.get_mapping_info().mapping,
+ matrix_free.get_dof_handler(),
+ Coordinates<dim>(),
+ src);
+
+ // 2) Construct remote evaluator...
+ FERemoteEvaluation<dim, dim, Number> phi_remote_eval(
+ remote_communicator, matrix_free.get_dof_handler());
+
+ // ...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) {
+ FEFacePointEvaluation<dim, dim, dim, Number> phi(
+ nm_mapping_info, matrix_free.get_dof_handler().get_fe());
+ auto phi_remote = phi_remote_eval.get_data_accessor();
+
+
+ for (unsigned int face = face_range.first; face < face_range.second;
+ ++face)
+ {
+ if (matrix_free.get_boundary_id(face) == 1)
+ {
+ for (unsigned int v = 0;
+ v < matrix_free.n_active_entries_per_face_batch(face);
+ ++v)
+ {
+ constexpr unsigned int n_lanes =
+ VectorizedArray<Number>::size();
+ phi.reinit(face * n_lanes + v);
+ phi_remote.reinit(face * n_lanes + v);
+
+ for (unsigned int q : phi.quadrature_point_indices())
+ {
+ const auto point = phi.real_point(q);
+ Tensor<1, dim, Number> temp;
+ for (unsigned int i = 0; i < dim; ++i)
+ temp[i] = point[i];
+
+ 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<VectorType, VectorType>(
+ {}, {}, boundary_function, dst, src, true);
+}
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+ MPILogInitAll all;
+
+ // Construct non-matching triangulation.
+ const double length = 1.0;
+ Triangulation<dim> tria_left;
+ const unsigned int subdiv_left = 11;
+ GridGenerator::subdivided_hyper_rectangle(tria_left,
+ {subdiv_left, 2 * subdiv_left},
+ {0.0, 0.0},
+ {0.525 * length, length});
+
+ for (const auto &face : tria_left.active_face_iterators())
+ if (face->at_boundary())
+ {
+ face->set_boundary_id(0);
+ if (face->center()[0] > 0.525 * length - 1e-6)
+ face->set_boundary_id(1);
+ }
+
+ Triangulation<dim> tria_right;
+ const unsigned int subdiv_right = 4;
+ GridGenerator::subdivided_hyper_rectangle(tria_right,
+ {subdiv_right, 2 * subdiv_right},
+ {0.525 * length, 0.0},
+ {length, length});
+
+ for (const auto &face : tria_right.active_face_iterators())
+ if (face->at_boundary())
+ face->set_boundary_id(0);
+
+ // Merge triangulations with tolerance 0 to ensure no vertices
+ // are merged.
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::merge_triangulations(tria_left,
+ tria_right,
+ tria,
+ /*tolerance*/ 0.,
+ /*copy_manifold_ids*/ false,
+ /*copy_boundary_ids*/ true);
+
+ // Create DoFHandler
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(3), dim));
+
+ // Constraints are not considered in this test
+ AffineConstraints<Number> constraints;
+ constraints.close();
+
+ // Setup MatrixFree
+ typename MatrixFree<dim, Number>::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;
+
+
+ MatrixFree<dim, Number> matrix_free;
+ matrix_free.reinit(
+ MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(4), data);
+
+ // 2) test point-to-point
+ test_point_to_point(matrix_free);
+
+ deallog << "OK" << std::endl;
+
+ // 3) test Nitsche-type
+ test_nitsche_type_mortaring(matrix_free);
+
+ deallog << "OK" << std::endl;
+
+ return 0;
+}
projects / dealii.git / commitdiff