16#ifndef dealii_matrix_free_fe_remote_evaluation_h
17#define dealii_matrix_free_fe_remote_evaluation_h
41 template <
int dim,
int n_components,
typename value_type_>
44 using value_type =
typename internal::FEPointEvaluation::
45 EvaluatorTypeTraits<dim, dim, n_components, value_type_>::value_type;
48 EvaluatorTypeTraits<dim, dim, n_components, value_type_>::
83 const unsigned int face_number)
const;
104 std::vector<unsigned int>
ptrs;
112 template <
int dim,
int n_components,
typename value_type_>
166 reinit(
const unsigned int index_0,
const unsigned int index_1);
205 std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>>
rpe;
218 std::vector<std::pair<unsigned int, unsigned int>>
234 std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>>
rpe;
246 std::vector<unsigned int>
262 std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>>
rpe;
269 std::pair<typename Triangulation<dim>::cell_iterator,
unsigned int>>
278 std::pair<typename Triangulation<dim>::cell_iterator,
unsigned int>>
297 const std::pair<unsigned int, unsigned int> &face_batch_range,
308 const std::pair<unsigned int, unsigned int> &face_range,
317 template <
typename Iterator>
329 template <
int n_components,
330 typename PrecomputedEvaluationDataType,
335 PrecomputedEvaluationDataType &dst,
336 const MeshType &mesh,
337 const VectorType &src,
339 const unsigned int first_selected_component,
358 std::vector<std::variant<FERemoteCommunicationObjectEntityBatches<dim>,
374 template <
typename T1,
typename T2>
378 const std::vector<T2> &src,
379 const std::vector<unsigned int> &indices);
385 template <
typename T1,
typename T2>
390 const std::vector<T2> &src,
392 unsigned int>> &cell_face_nos);
398 template <
typename T1,
typename T2>
402 const std::vector<T2> &src,
403 const std::vector<std::pair<unsigned int, unsigned int>>
404 &batch_id_n_entities);
410 template <
typename T1, std::
size_t n_lanes>
413 const unsigned int v,
419 template <
typename T1,
int rank_, std::
size_t n_lanes,
int dim_>
422 const unsigned int v,
428 template <
typename T1,
438 const unsigned int v,
445 template <
typename T1,
typename T2>
484 &non_matching_faces_marked_vertices,
485 const unsigned int quad_no = 0,
486 const unsigned int dof_no = 0,
487 const double tolerance = 1e-9);
517 &non_matching_faces_marked_vertices,
518 const unsigned int n_q_pnts_1D,
519 const unsigned int dof_no = 0,
521 const double tolerance = 1e-9);
537template <
int dim,
int n_components,
typename value_type>
554 template <
typename MeshType>
556 const MeshType &mesh,
569 template <
typename VectorType>
654 const unsigned int face_number)
const
668 return ptrs[face_index];
678 template <
int dim,
int n_components,
typename value_type_>
685 , data_offset(
numbers::invalid_unsigned_int)
688 template <
int dim,
int n_components,
typename value_type_>
691 value_type_>::value_type
693 const unsigned int q)
const
698 return data.values[data_offset + q];
701 template <
int dim,
int n_components,
typename value_type_>
710 return data.gradients[data_offset + q];
713 template <
int dim,
int n_components,
typename value_type_>
716 const unsigned int index)
718 data_offset = view.get_shift(
index);
721 template <
int dim,
int n_components,
typename value_type_>
724 const unsigned int index_0,
725 const unsigned int index_1)
727 data_offset = view.get_shift(index_0, index_1);
735std::vector<std::pair<unsigned int, unsigned int>>
739 return batch_id_n_entities;
743std::vector<unsigned int>
750std::vector<std::pair<typename Triangulation<dim>::cell_iterator,
unsigned int>>
753 return cell_face_nos;
763 const std::pair<unsigned int, unsigned int> &face_batch_range,
767 communication_objects.clear();
768 for (
const auto &co : comm_objects)
769 communication_objects.push_back(co);
772 const unsigned int n_cells = quadrature_vector.size();
773 AssertDimension(n_cells, face_batch_range.second - face_batch_range.first);
776 view.start = face_batch_range.first;
778 view.ptrs.resize(n_cells + 1);
781 for (
unsigned int face = 0; face < n_cells; ++face)
783 view.ptrs[face + 1] = view.ptrs[face] + quadrature_vector[face].size();
791 const std::pair<unsigned int, unsigned int> &face_range,
795 communication_objects.clear();
796 for (
const auto &co : comm_objects)
797 communication_objects.push_back(co);
799 const unsigned int n_faces = quadrature_vector.size();
803 view.start = face_range.first;
805 view.ptrs.resize(n_faces + 1);
808 for (
unsigned int face = 0; face < n_faces; ++face)
809 view.ptrs[face + 1] = view.ptrs[face] + quadrature_vector[face].size();
813template <
typename Iterator>
821 communication_objects.clear();
822 for (
const auto &co : comm_objects)
823 communication_objects.push_back(co);
825 const unsigned int n_cells = quadrature_vector.size();
827 std::distance(cell_iterator_range.
begin(),
828 cell_iterator_range.
end()));
831 auto &cell_ptrs = view.ptrs_ptrs;
832 auto &face_ptrs = view.ptrs;
835 cell_ptrs.resize(n_cells);
836 unsigned int n_faces = 0;
837 for (
const auto &cell : cell_iterator_range)
839 cell_ptrs[cell->active_cell_index()] = n_faces;
840 n_faces += cell->n_faces();
843 face_ptrs.resize(n_faces + 1);
845 for (
const auto &cell : cell_iterator_range)
847 for (
const auto &f : cell->face_indices())
849 const unsigned int face_index =
850 cell_ptrs[cell->active_cell_index()] + f;
852 face_ptrs[face_index + 1] =
853 face_ptrs[face_index] +
854 quadrature_vector[cell->active_cell_index()][f].size();
860template <
int n_components,
861 typename PrecomputedEvaluationDataType,
866 PrecomputedEvaluationDataType &dst,
867 const MeshType &mesh,
868 const VectorType &src,
870 const unsigned int first_selected_component,
873 const bool has_ghost_elements = src.has_ghost_elements();
875 if (has_ghost_elements ==
false)
876 src.update_ghost_values();
879 for (
const auto &communication_object : communication_objects)
884 [&](
const auto &obj) {
885 CopyInstructions::copy_data(
888 VectorTools::point_values<n_components>(
889 *obj.rpe, mesh, src, vec_flags, first_selected_component),
890 obj.get_communication_object_pntrs());
892 communication_object);
898 [&](
const auto &obj) {
899 CopyInstructions::copy_data(
902 VectorTools::point_gradients<n_components>(
903 *obj.rpe, mesh, src, vec_flags, first_selected_component),
904 obj.get_communication_object_pntrs());
906 communication_object);
912 if (has_ghost_elements ==
false)
913 src.zero_out_ghost_values();
924template <
typename T1,
typename T2>
929 const std::vector<T2> &src,
930 const std::vector<unsigned int> &indices)
935 for (
const auto idx : indices)
948template <
typename T1,
typename T2>
953 const std::vector<T2> &src,
955 unsigned int>> &cell_face_nos)
960 for (
const auto &[cell, f] : cell_face_nos)
962 for (
unsigned int j =
view.
get_shift(cell->active_cell_index(), f);
975template <
typename T1,
typename T2>
980 const std::vector<T2> &src,
981 const std::vector<std::pair<unsigned int, unsigned int>> &batch_id_n_entities)
986 for (
const auto &[batch_id, n_entries] : batch_id_n_entities)
988 for (
unsigned int v = 0; v < n_entries; ++v)
996 copy_data_entries(dst[j], v, src[c]);
1002template <
typename T1, std::
size_t n_lanes>
1006 const unsigned int v,
1015template <
typename T1,
int rank_, std::
size_t n_lanes,
int dim_>
1019 const unsigned int v,
1024 if constexpr (rank_ == 1)
1026 for (
unsigned int i = 0; i < dim_; ++i)
1031 for (
unsigned int i = 0; i < rank_; ++i)
1032 for (
unsigned int j = 0; j < dim_; ++j)
1033 dst[i][j][v] = src[i][j];
1038template <
typename T1,
1040 std::size_t n_lanes,
1048 const unsigned int v,
1051 if constexpr (rank_ == 1)
1053 for (
unsigned int i = 0; i < n_components_; ++i)
1054 copy_data(dst[i], v, src[i]);
1058 for (
unsigned int i = 0; i < rank_; ++i)
1059 for (
unsigned int j = 0; j < n_components_; ++j)
1060 dst[i][j][v] = src[i][j];
1065template <
typename T1,
typename T2>
1074 "copy_data_entries() not implemented for given arguments."));
1081 template <
int dim,
typename Number,
typename VectorizedArrayType>
1087 &non_matching_faces_marked_vertices,
1088 const unsigned int quad_no,
1089 const unsigned int dof_no,
1090 const double tolerance)
1109 std::vector<FERemoteCommunicationObjectEntityBatches<dim>> comm_objects;
1116 std::vector<Quadrature<dim>> global_quadrature_vector(
1121 const auto face_batch_range =
1127 for (
const auto &[nm_face, marked_vertices] :
1128 non_matching_faces_marked_vertices)
1133 auto rpe = std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(
1134 tolerance,
false, 0, marked_vertices);
1137 std::vector<std::pair<unsigned int, unsigned int>>
1138 face_batch_id_n_faces;
1141 std::vector<Point<dim>> points;
1151 for (
unsigned int bface = 0;
1152 bface < face_batch_range.second - face_batch_range.first;
1155 const unsigned int face = face_batch_range.first + bface;
1164 const unsigned int n_faces =
1166 face_batch_id_n_faces.emplace_back(
1167 std::make_pair(face, n_faces));
1171 for (
unsigned int v = 0; v < n_faces; ++v)
1173 for (
unsigned int q : phi.quadrature_point_indices())
1175 const auto point = phi.quadrature_point(q);
1177 for (
unsigned int i = 0; i < dim; ++i)
1178 temp[i] = point[i][v];
1180 points.push_back(temp);
1187 Assert(global_quadrature_vector[bface].size() == 0,
1189 "Quadrature for given face already provided."));
1191 global_quadrature_vector[bface] =
1197 rpe->reinit(points, tria, mapping);
1198 Assert(rpe->all_points_found(),
1205 comm_objects.push_back(co);
1214 global_quadrature_vector);
1216 return remote_communicator;
1221 template <
int dim,
typename Number,
typename VectorizedArrayType>
1227 &non_matching_faces_marked_vertices,
1228 const unsigned int n_q_pnts_1D,
1229 const unsigned int dof_no,
1231 const double tolerance)
1239 std::pair<unsigned int, unsigned int> face_range =
1244 std::vector<
Quadrature<dim - 1>> global_quadrature_vector(
1256 std::vector<FERemoteCommunicationObject<dim>> comm_objects;
1260 std::vector<BoundingBox<dim>> local_boxes;
1262 if (cell->is_locally_owned())
1263 local_boxes.emplace_back(mapping.get_bounding_box(cell));
1266 const auto local_tree =
pack_rtree(local_boxes);
1269 std::vector<std::vector<BoundingBox<dim>>> global_bboxes(1);
1275 for (
const auto &[nm_face, marked_vertices] :
1276 non_matching_faces_marked_vertices)
1280 std::pair<typename Triangulation<dim>::cell_iterator,
unsigned int>>
1283 std::vector<unsigned int> indices;
1285 for (
unsigned int face = face_range.first; face < face_range.second;
1290 for (
unsigned int v = 0;
1296 cell_face_pairs.emplace_back(std::make_pair(c, f));
1297 indices.push_back(face * n_lanes + v);
1310 std::vector<std::vector<Point<dim>>> intersection_requests;
1311 for (
const auto &[cell, f] : cell_face_pairs)
1313 std::vector<Point<dim>>
vertices(cell->face(f)->n_vertices());
1314 std::copy_n(mapping.get_vertices(cell, f).begin(),
1315 cell->face(f)->n_vertices(),
1317 intersection_requests.emplace_back(
vertices);
1322 auto intersection_data =
1325 intersection_requests,
1331 std::vector<Quadrature<dim>> mapped_quadratures_recv_comp;
1334 std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>();
1337 .
template convert_to_distributed_compute_point_locations_internal<
1338 dim>(n_q_pnts_1D, tria, mapping, &mapped_quadratures_recv_comp),
1343 for (
unsigned int i = 0; i < intersection_requests.size(); ++i)
1345 const auto idx = indices[i];
1350 const auto &cell = std::get<0>(cell_face_pairs[i]);
1351 const auto &f = std::get<1>(cell_face_pairs[i]);
1353 std::vector<
Point<dim - 1>> q_points;
1354 std::vector<double> weights;
1355 for (
unsigned int ptr = intersection_data.recv_ptrs[i];
1356 ptr < intersection_data.recv_ptrs[i + 1];
1359 const auto &quad = mapped_quadratures_recv_comp[ptr];
1361 const auto &ps = quad.get_points();
1365 std::back_inserter(q_points),
1367 return mapping.project_real_point_to_unit_point_on_face(
1371 const auto &ws = quad.get_weights();
1372 weights.insert(weights.end(), ws.begin(), ws.end());
1376 Assert(global_quadrature_vector[idx].size() == 0,
1377 ExcMessage(
"Quadrature for given face already provided."));
1379 global_quadrature_vector[idx] = quad;
1386 comm_objects.push_back(co);
1394 std::make_pair(0, global_quadrature_vector.size()),
1395 global_quadrature_vector);
1397 if (nm_mapping_info !=
nullptr)
1400 std::pair<typename DoFHandler<dim>::cell_iterator,
unsigned int>>
1401 vector_face_accessors;
1407 unsigned int face_batch = 0;
1410 for (
unsigned int v = 0; v < n_lanes; ++v)
1413 vector_face_accessors.push_back(
1416 vector_face_accessors.push_back(
1425 for (
unsigned int v = 0; v < n_lanes; ++v)
1428 vector_face_accessors.push_back(
1431 vector_face_accessors.push_back(
1437 global_quadrature_vector);
1440 return remote_communicator;
1446template <
int dim,
int n_components,
typename value_type>
1447template <
typename MeshType>
1450 const MeshType &mesh,
1451 const unsigned int first_selected_component,
1454 , first_selected_component(first_selected_component)
1455 , evaluation_flags(evaluation_flags)
1460template <
int dim,
int n_components,
typename value_type>
1461template <
typename VectorType>
1464 const VectorType &src,
1470 comm->template update_ghost_values<n_components>(this->data,
1474 first_selected_component,
1477 else if (dof_handler)
1479 comm->template update_ghost_values<n_components>(this->data,
1483 first_selected_component,
1490template <
int dim,
int n_components,
typename value_type>
1496 return data_accessor;
1499template <
int dim,
int n_components,
typename value_type>
1507template <
int dim,
int n_components,
typename value_type>
1512 this->dof_handler = &dof_handler;
void resize(const size_type new_size)
static void copy_data(const internal::PrecomputedEvaluationDataView &view, AlignedVector< T1 > &dst, const std::vector< T2 > &src, const std::vector< unsigned int > &indices)
static void copy_data_entries(VectorizedArray< T1, n_lanes > &dst, const unsigned int v, const T1 &src)
void reinit_faces(const std::vector< FERemoteCommunicationObjectEntityBatches< dim > > &comm_objects, const std::pair< unsigned int, unsigned int > &face_batch_range, const std::vector< Quadrature< dim > > &quadrature_vector)
const internal::PrecomputedEvaluationDataView & get_view() const
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
internal::PrecomputedEvaluationDataView view
std::vector< std::variant< FERemoteCommunicationObjectEntityBatches< dim >, FERemoteCommunicationObject< dim >, FERemoteCommunicationObjectTwoLevel< dim > > > communication_objects
SmartPointer< const DoFHandler< dim > > dof_handler
FERemoteEvaluation(const FERemoteEvaluationCommunicator< dim > &comm, const MeshType &mesh, const unsigned int first_selected_component=0, const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags=VectorTools::EvaluationFlags::avg)
SmartPointer< const Triangulation< dim > > tria
void set_mesh(const Triangulation< dim > &tria)
SmartPointer< const FERemoteEvaluationCommunicator< dim > > comm
internal::PrecomputedEvaluationData< dim, n_components, value_type > data
const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags
const unsigned int first_selected_component
internal::PrecomputedEvaluationDataAccessor< dim, n_components, value_type > get_data_accessor() const
void gather_evaluate(const VectorType &src, const EvaluationFlags::EvaluationFlags flags)
IteratorOverIterators end() const
IteratorOverIterators begin()
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
unsigned int n_inner_face_batches() const
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
unsigned int n_boundary_face_batches() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
void reinit_faces(const ContainerType &cell_iterator_range, const std::vector< std::vector< Quadrature< dim - 1 > > > &quadrature_vector, const unsigned int n_unfiltered_cells=numbers::invalid_unsigned_int)
Triangulation< dim, spacedim > & get_triangulation()
static constexpr std::size_t size()
const gradient_type get_gradient(const unsigned int q) const
PrecomputedEvaluationDataAccessor(const PrecomputedEvaluationData< dim, n_components, value_type_ > &data, const PrecomputedEvaluationDataView &view)
const PrecomputedEvaluationDataView & view
typename PrecomputedEvaluationData< dim, n_components, value_type_ >::value_type value_type
const PrecomputedEvaluationData< dim, n_components, value_type_ > & data
const value_type get_value(const unsigned int q) const
void reinit(const unsigned int index)
typename PrecomputedEvaluationData< dim, n_components, value_type_ >::gradient_type gradient_type
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NOT_IMPLEMENTED()
EvaluationFlags
The EvaluationFlags enum.
FERemoteEvaluationCommunicator< dim > compute_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)
FERemoteEvaluationCommunicator< dim > compute_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)
static const unsigned int invalid_unsigned_int
*braid_SplitCommworld & comm
std::vector< BoundingBox< boost::geometry::dimension< typename Rtree::indexable_type >::value > > extract_rtree_level(const Rtree &tree, const unsigned int level)
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
std::vector< std::pair< unsigned int, unsigned int > > batch_id_n_entities
std::vector< std::pair< unsigned int, unsigned int > > get_communication_object_pntrs() const
std::shared_ptr< Utilities::MPI::RemotePointEvaluation< dim > > rpe
std::shared_ptr< Utilities::MPI::RemotePointEvaluation< dim > > rpe
std::vector< std::pair< typename Triangulation< dim >::cell_iterator, unsigned int > > cell_face_nos
std::vector< std::pair< typename Triangulation< dim >::cell_iterator, unsigned int > > get_communication_object_pntrs() const
std::shared_ptr< Utilities::MPI::RemotePointEvaluation< dim > > rpe
std::vector< unsigned int > indices
std::vector< unsigned int > get_communication_object_pntrs() const
unsigned int get_shift(const unsigned int index) const
unsigned int size() const
std::vector< unsigned int > ptrs_ptrs
std::vector< unsigned int > ptrs
AlignedVector< value_type > values
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, dim, n_components, value_type_ >::value_type value_type
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, dim, n_components, value_type_ >::real_gradient_type gradient_type
AlignedVector< gradient_type > gradients