Reference documentation for deal.II version GIT relicensing-136-gb80d0be4af 2024-03-18 08:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
quadrature_point_data.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_quadrature_point_data_h
16#define dealii_quadrature_point_data_h
17
18#include <deal.II/base/config.h>
19
22
24
25#include <deal.II/fe/fe.h>
26#include <deal.II/fe/fe_tools.h>
27
28#include <deal.II/grid/tria.h>
30
31#include <deal.II/lac/vector.h>
32
33#include <map>
34#include <optional>
35#include <type_traits>
36#include <vector>
37
39
62template <typename CellIteratorType, typename DataType>
64{
65public:
69 CellDataStorage() = default;
70
74 ~CellDataStorage() override = default;
75
101 template <typename T = DataType>
102 void
103 initialize(const CellIteratorType &cell,
104 const unsigned int number_of_data_points_per_cell);
105
111 template <typename T = DataType>
112 void
113 initialize(const CellIteratorType &cell_start,
114 const CellIteratorType &cell_end,
115 const unsigned int number_of_data_points_per_cell);
116
126 bool
127 erase(const CellIteratorType &cell);
128
132 void
134
149 template <typename T = DataType>
150 std::vector<std::shared_ptr<T>>
151 get_data(const CellIteratorType &cell);
152
167 template <typename T = DataType>
168 std::vector<std::shared_ptr<const T>>
169 get_data(const CellIteratorType &cell) const;
170
187 template <typename T = DataType>
188 std::optional<std::vector<std::shared_ptr<T>>>
189 try_get_data(const CellIteratorType &cell);
190
207 template <typename T = DataType>
208 std::optional<std::vector<std::shared_ptr<const T>>>
209 try_get_data(const CellIteratorType &cell) const;
210
211private:
215 static constexpr unsigned int dimension =
216 CellIteratorType::AccessorType::dimension;
217
221 static constexpr unsigned int space_dimension =
222 CellIteratorType::AccessorType::space_dimension;
223
232
238 std::map<CellId, std::vector<std::shared_ptr<DataType>>> map;
239
245 "Cell data is being retrieved with a type which is different than the type used to initialize it");
246
252 "The provided cell iterator does not belong to the triangulation that corresponds to the CellDataStorage object.");
253};
254
255
277{
278public:
283
288
294 virtual unsigned int
295 number_of_values() const = 0;
296
307 virtual void
308 pack_values(std::vector<double> &values) const = 0;
309
318 virtual void
319 unpack_values(const std::vector<double> &values) = 0;
320};
321
322
323#ifdef DEAL_II_WITH_P4EST
324namespace parallel
325{
326 namespace distributed
327 {
437 template <int dim, typename DataType>
439 {
440 public:
441 static_assert(
442 std::is_base_of_v<TransferableQuadraturePointData, DataType>,
443 "User's DataType class should be derived from TransferableQuadraturePointData");
444
450
462 const Quadrature<dim> &mass_quadrature,
463 const Quadrature<dim> &data_quadrature);
464
475 void
479
490 void
492
493 private:
499 std::vector<char>
502 &cell,
503 const CellStatus status);
504
510 void
513 &cell,
514 const CellStatus status,
515 const boost::iterator_range<std::vector<char>::const_iterator>
516 &data_range);
517
521 const std::unique_ptr<const FiniteElement<dim>> projection_fe;
522
528
532 const unsigned int n_q_points;
533
539
545
553
558
564
569 unsigned int handle;
570
575
581 };
582
583 } // namespace distributed
584
585} // namespace parallel
586
587#endif
588
591#ifndef DOXYGEN
592
593// ------------------- inline and template functions ----------------
594
595//--------------------------------------------------------------------
596// CellDataStorage
597//--------------------------------------------------------------------
598
599template <typename CellIteratorType, typename DataType>
600template <typename T>
601inline void
603 const CellIteratorType &cell,
604 const unsigned int n_q_points)
605{
606 static_assert(std::is_base_of_v<DataType, T>,
607 "User's T class should be derived from user's DataType class");
608 // The first time this method is called, it has to initialize the reference
609 // to the triangulation object
610 if (!tria)
611 tria = &cell->get_triangulation();
612 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
613
614 const auto key = cell->id();
615 if (map.find(key) == map.end())
616 {
617 map[key] = std::vector<std::shared_ptr<DataType>>(n_q_points);
618 // we need to initialize one-by-one as the std::vector<>(q, T())
619 // will end with a single same T object stored in each element of the
620 // vector:
621 const auto it = map.find(key);
622 for (unsigned int q = 0; q < n_q_points; ++q)
623 it->second[q] = std::make_shared<T>();
624 }
625}
626
627
628
629template <typename CellIteratorType, typename DataType>
630template <typename T>
631inline void
633 const CellIteratorType &cell_start,
634 const CellIteratorType &cell_end,
635 const unsigned int number)
636{
637 for (CellIteratorType it = cell_start; it != cell_end; ++it)
638 if (it->is_locally_owned())
639 initialize<T>(it, number);
640}
641
642
643
644template <typename CellIteratorType, typename DataType>
645inline bool
646CellDataStorage<CellIteratorType, DataType>::erase(const CellIteratorType &cell)
647{
648 const auto key = cell->id();
649 const auto it = map.find(key);
650 if (it == map.end())
651 return false;
652 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
653 for (unsigned int i = 0; i < it->second.size(); ++i)
654 {
655 Assert(
656 it->second[i].unique(),
658 "Can not erase the cell data multiple objects reference its data."));
659 }
660
661 return (map.erase(key) == 1);
662}
663
664
665
666template <typename CellIteratorType, typename DataType>
667inline void
669{
670 // Do not call
671 // map.clear();
672 // as we want to be sure no one uses the stored objects. Loop manually:
673 auto it = map.begin();
674 while (it != map.end())
675 {
676 // loop over all objects and see if no one is using them
677 for (unsigned int i = 0; i < it->second.size(); ++i)
678 {
679 Assert(
680 it->second[i].unique(),
682 "Can not erase the cell data, multiple objects reference it."));
683 }
684 it = map.erase(it);
685 }
686}
687
688
689
690template <typename CellIteratorType, typename DataType>
691template <typename T>
692inline std::vector<std::shared_ptr<T>>
694 const CellIteratorType &cell)
695{
696 static_assert(std::is_base_of_v<DataType, T>,
697 "User's T class should be derived from user's DataType class");
698 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
699
700 const auto it = map.find(cell->id());
701 Assert(it != map.end(), ExcMessage("Could not find data for the cell"));
702
703 // It would be nice to have a specialized version of this function for
704 // T==DataType. However explicit (i.e full) specialization of a member
705 // template is only allowed when the enclosing class is also explicitly (i.e
706 // fully) specialized. Thus, stick with copying of shared pointers even when
707 // the T==DataType:
708 std::vector<std::shared_ptr<T>> res(it->second.size());
709 for (unsigned int q = 0; q < res.size(); ++q)
710 {
711 res[q] = std::dynamic_pointer_cast<T>(it->second[q]);
712 Assert(res[q], ExcCellDataTypeMismatch());
713 }
714 return res;
715}
716
717
718
719template <typename CellIteratorType, typename DataType>
720template <typename T>
721inline std::vector<std::shared_ptr<const T>>
723 const CellIteratorType &cell) const
724{
725 static_assert(std::is_base_of_v<DataType, T>,
726 "User's T class should be derived from user's DataType class");
727 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
728
729 const auto it = map.find(cell->id());
730 Assert(it != map.end(), ExcMessage("Could not find QP data for the cell"));
731
732 // Cast base class to the desired class. This has to be done irrespectively of
733 // T==DataType as we need to return shared_ptr<const T> to make sure the user
734 // does not modify the content of QP objects
735 std::vector<std::shared_ptr<const T>> res(it->second.size());
736 for (unsigned int q = 0; q < res.size(); ++q)
737 {
738 res[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
739 Assert(res[q], ExcCellDataTypeMismatch());
740 }
741 return res;
742}
743
744template <typename CellIteratorType, typename DataType>
745template <typename T>
746inline std::optional<std::vector<std::shared_ptr<T>>>
748 const CellIteratorType &cell)
749{
750 static_assert(std::is_base_of_v<DataType, T>,
751 "User's T class should be derived from user's DataType class");
752 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
753
754 const auto it = map.find(cell->id());
755 if (it != map.end())
756 {
757 // Cast base class to the desired class. This has to be done
758 // irrespectively of T==DataType as we need to return
759 // shared_ptr<const T> to make sure the user
760 // does not modify the content of QP objects
761 std::vector<std::shared_ptr<T>> result(it->second.size());
762 for (unsigned int q = 0; q < result.size(); ++q)
763 {
764 result[q] = std::dynamic_pointer_cast<T>(it->second[q]);
765 Assert(result[q], ExcCellDataTypeMismatch());
766 }
767 return {result};
768 }
769 else
770 {
771 return {};
772 }
773}
774
775template <typename CellIteratorType, typename DataType>
776template <typename T>
777inline std::optional<std::vector<std::shared_ptr<const T>>>
779 const CellIteratorType &cell) const
780{
781 static_assert(std::is_base_of_v<DataType, T>,
782 "User's T class should be derived from user's DataType class");
783 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
784
785 const auto it = map.find(cell->id());
786 if (it != map.end())
787 {
788 // Cast base class to the desired class. This has to be done
789 // irrespectively of T==DataType as we need to return
790 // shared_ptr<const T> to make sure the user
791 // does not modify the content of QP objects
792 std::vector<std::shared_ptr<const T>> result(it->second.size());
793 for (unsigned int q = 0; q < result.size(); ++q)
794 {
795 result[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
796 Assert(result[q], ExcCellDataTypeMismatch());
797 }
798 return {result};
799 }
800 else
801 {
802 return {};
803 }
804}
805
806//--------------------------------------------------------------------
807// ContinuousQuadratureDataTransfer
808//--------------------------------------------------------------------
809
810
811/*
812 * Pack cell data of type @p DataType stored using @p data_storage in @p cell
813 * at each quadrature point to @p matrix_data. Here @p matrix_data is a matrix
814 * whose first index corresponds to different quadrature points on the cell
815 * whereas the second index represents different values stored at each
816 * quadrature point in the DataType class.
817 */
818template <typename CellIteratorType, typename DataType>
819inline void
820pack_cell_data(const CellIteratorType &cell,
822 FullMatrix<double> &matrix_data)
823{
824 static_assert(std::is_base_of_v<TransferableQuadraturePointData, DataType>,
825 "User's DataType class should be derived from QPData");
826
827 if (const auto qpd = data_storage->try_get_data(cell))
828 {
829 const unsigned int m = qpd->size();
830 Assert(m > 0, ExcInternalError());
831 const unsigned int n = (*qpd)[0]->number_of_values();
832 matrix_data.reinit(m, n);
833
834 std::vector<double> single_qp_data(n);
835 for (unsigned int q = 0; q < m; ++q)
836 {
837 (*qpd)[q]->pack_values(single_qp_data);
838 AssertDimension(single_qp_data.size(), n);
839
840 for (unsigned int i = 0; i < n; ++i)
841 matrix_data(q, i) = single_qp_data[i];
842 }
843 }
844 else
845 {
846 matrix_data.reinit({0, 0});
847 }
848}
849
850
851
852/*
853 * the opposite of the pack function above.
854 */
855template <typename CellIteratorType, typename DataType>
856inline void
857unpack_to_cell_data(const CellIteratorType &cell,
858 const FullMatrix<double> &values_at_qp,
860{
861 static_assert(std::is_base_of_v<TransferableQuadraturePointData, DataType>,
862 "User's DataType class should be derived from QPData");
863
864 if (const auto qpd = data_storage->try_get_data(cell))
865 {
866 const unsigned int n = values_at_qp.n();
867 AssertDimension((*qpd)[0]->number_of_values(), n);
868
869 std::vector<double> single_qp_data(n);
870 AssertDimension(qpd->size(), values_at_qp.m());
871
872 for (unsigned int q = 0; q < qpd->size(); ++q)
873 {
874 for (unsigned int i = 0; i < n; ++i)
875 single_qp_data[i] = values_at_qp(q, i);
876 (*qpd)[q]->unpack_values(single_qp_data);
877 }
878 }
879}
880
881
882# ifdef DEAL_II_WITH_P4EST
883
884namespace parallel
885{
886 namespace distributed
887 {
888 template <int dim, typename DataType>
891 const Quadrature<dim> &lhs_quadrature,
892 const Quadrature<dim> &rhs_quadrature)
893 : projection_fe(
894 std::unique_ptr<const FiniteElement<dim>>(projection_fe_.clone()))
895 , data_size_in_bytes(0)
896 , n_q_points(rhs_quadrature.size())
897 , project_to_fe_matrix(projection_fe->n_dofs_per_cell(), n_q_points)
898 , project_to_qp_matrix(n_q_points, projection_fe->n_dofs_per_cell())
899 , handle(numbers::invalid_unsigned_int)
900 , data_storage(nullptr)
901 , triangulation(nullptr)
902 {
903 Assert(
904 projection_fe->n_components() == 1,
906 "ContinuousQuadratureDataTransfer requires scalar FiniteElement"));
907
909 *projection_fe.get(),
910 lhs_quadrature,
911 rhs_quadrature,
913
915 *projection_fe.get(), rhs_quadrature, project_to_qp_matrix);
916 }
917
918
919
920 template <int dim, typename DataType>
921 inline void
926 {
927 Assert(data_storage == nullptr,
928 ExcMessage("This function can be called only once"));
929 triangulation = &tr_;
930 data_storage = &data_storage_;
931
932 handle = triangulation->register_data_attach(
933 [this](const typename parallel::distributed::Triangulation<
934 dim>::cell_iterator &cell,
935 const CellStatus status) {
936 return this->pack_function(cell, status);
937 },
938 /*returns_variable_size_data=*/true);
939 }
940
941
942
943 template <int dim, typename DataType>
944 inline void
946 {
947 triangulation->notify_ready_to_unpack(
948 handle,
949 [this](const typename parallel::distributed::Triangulation<
950 dim>::cell_iterator &cell,
951 const CellStatus status,
952 const boost::iterator_range<std::vector<char>::const_iterator>
953 &data_range) {
954 this->unpack_function(cell, status, data_range);
955 });
956
957 // invalidate the pointers
958 data_storage = nullptr;
959 triangulation = nullptr;
960 }
961
962
963
964 template <int dim, typename DataType>
965 inline std::vector<char>
968 &cell,
969 const CellStatus /*status*/)
970 {
971 pack_cell_data(cell, data_storage, matrix_quadrature);
972
973 // project to FE
974 const unsigned int number_of_values = matrix_quadrature.n();
975 matrix_dofs.reinit(project_to_fe_matrix.m(), number_of_values);
976 if (number_of_values > 0)
977 project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
978
979 return Utilities::pack(matrix_dofs, /*allow_compression=*/false);
980 }
981
982
983
984 template <int dim, typename DataType>
985 inline void
988 &cell,
989 const CellStatus status,
990 const boost::iterator_range<std::vector<char>::const_iterator>
991 &data_range)
992 {
995 (void)status;
996
997 matrix_dofs =
998 Utilities::unpack<FullMatrix<double>>(data_range.begin(),
999 data_range.end(),
1000 /*allow_compression=*/false);
1001 const unsigned int number_of_values = matrix_dofs.n();
1002 if (number_of_values == 0)
1003 return;
1004
1005 matrix_quadrature.reinit(n_q_points, number_of_values);
1006
1007 if (cell->has_children())
1008 {
1009 // we need to first use prolongation matrix to get dofvalues on child
1010 // cells based on dofvalues stored in the parent's data_store
1011 matrix_dofs_child.reinit(projection_fe->n_dofs_per_cell(),
1012 number_of_values);
1013 for (unsigned int child = 0; child < cell->n_children(); ++child)
1014 if (cell->child(child)->is_locally_owned())
1015 {
1016 projection_fe
1017 ->get_prolongation_matrix(child, cell->refinement_case())
1018 .mmult(matrix_dofs_child, matrix_dofs);
1019
1020 // now we do the usual business of evaluating FE on quadrature
1021 // points:
1022 project_to_qp_matrix.mmult(matrix_quadrature,
1023 matrix_dofs_child);
1024
1025 // finally, put back into the map:
1026 unpack_to_cell_data(cell->child(child),
1027 matrix_quadrature,
1028 data_storage);
1029 }
1030 }
1031 else
1032 {
1033 // if there are no children, evaluate FE field at
1034 // rhs_quadrature points.
1035 project_to_qp_matrix.mmult(matrix_quadrature, matrix_dofs);
1036
1037 // finally, put back into the map:
1038 unpack_to_cell_data(cell, matrix_quadrature, data_storage);
1039 }
1040 }
1041
1042 } // namespace distributed
1043
1044} // namespace parallel
1045
1046# endif // DEAL_II_WITH_P4EST
1047
1048#endif // DOXYGEN
1050
1051#endif
CellStatus
Definition cell_status.h:31
@ children_will_be_coarsened
std::optional< std::vector< std::shared_ptr< const T > > > try_get_data(const CellIteratorType &cell) const
CellDataStorage()=default
void initialize(const CellIteratorType &cell, const unsigned int number_of_data_points_per_cell)
bool erase(const CellIteratorType &cell)
static constexpr unsigned int space_dimension
SmartPointer< const Triangulation< dimension, space_dimension >, CellDataStorage< CellIteratorType, DataType > > tria
std::vector< std::shared_ptr< const T > > get_data(const CellIteratorType &cell) const
void initialize(const CellIteratorType &cell_start, const CellIteratorType &cell_end, const unsigned int number_of_data_points_per_cell)
std::optional< std::vector< std::shared_ptr< T > > > try_get_data(const CellIteratorType &cell)
~CellDataStorage() override=default
std::map< CellId, std::vector< std::shared_ptr< DataType > > > map
std::vector< std::shared_ptr< T > > get_data(const CellIteratorType &cell)
static constexpr unsigned int dimension
size_type n() const
size_type m() const
virtual ~TransferableQuadraturePointData()=default
virtual unsigned int number_of_values() const =0
virtual void unpack_values(const std::vector< double > &values)=0
virtual void pack_values(std::vector< double > &values) const =0
Triangulation< dim, spacedim > & get_triangulation()
const std::unique_ptr< const FiniteElement< dim > > projection_fe
ContinuousQuadratureDataTransfer(const FiniteElement< dim > &projection_fe, const Quadrature< dim > &mass_quadrature, const Quadrature< dim > &data_quadrature)
parallel::distributed::Triangulation< dim > * triangulation
typename parallel::distributed::Triangulation< dim >::cell_iterator CellIteratorType
std::vector< char > pack_function(const typename parallel::distributed::Triangulation< dim >::cell_iterator &cell, const CellStatus status)
void prepare_for_coarsening_and_refinement(parallel::distributed::Triangulation< dim > &tria, CellDataStorage< CellIteratorType, DataType > &data_storage)
void unpack_function(const typename parallel::distributed::Triangulation< dim >::cell_iterator &cell, const CellStatus status, const boost::iterator_range< std::vector< char >::const_iterator > &data_range)
CellDataStorage< CellIteratorType, DataType > * data_storage
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcCellDataTypeMismatch()
#define AssertDimension(dim1, dim2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcTriangulationMismatch()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition tria.h:269
void compute_interpolation_to_quadrature_points_matrix(const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, FullMatrix< double > &I_q)
void compute_projection_from_quadrature_points_matrix(const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &lhs_quadrature, const Quadrature< dim > &rhs_quadrature, FullMatrix< double > &X)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition utilities.h:1364
STL namespace.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation