deal.II version GIT relicensing-2202-g5ed228918a 2024-12-01 17:10:00+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 - 2024 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
114 const CellIteratorType &cell_start,
116 const unsigned int number_of_data_points_per_cell);
117
127 bool
128 erase(const CellIteratorType &cell);
129
133 void
135
150 template <typename T = DataType>
151 std::vector<std::shared_ptr<T>>
152 get_data(const CellIteratorType &cell);
153
168 template <typename T = DataType>
169 std::vector<std::shared_ptr<const T>>
170 get_data(const CellIteratorType &cell) const;
171
188 template <typename T = DataType>
189 std::optional<std::vector<std::shared_ptr<T>>>
190 try_get_data(const CellIteratorType &cell);
191
208 template <typename T = DataType>
209 std::optional<std::vector<std::shared_ptr<const T>>>
210 try_get_data(const CellIteratorType &cell) const;
211
212private:
216 static constexpr unsigned int dimension =
217 CellIteratorType::AccessorType::dimension;
218
222 static constexpr unsigned int space_dimension =
223 CellIteratorType::AccessorType::space_dimension;
224
233
239 std::map<CellId, std::vector<std::shared_ptr<DataType>>> map;
240
246 "Cell data is being retrieved with a type which is different than the type used to initialize it");
247
253 "The provided cell iterator does not belong to the triangulation that corresponds to the CellDataStorage object.");
254};
255
256
278{
279public:
284
289
295 virtual unsigned int
296 number_of_values() const = 0;
297
308 virtual void
309 pack_values(std::vector<double> &values) const = 0;
310
319 virtual void
320 unpack_values(const std::vector<double> &values) = 0;
321};
322
323
324#ifdef DEAL_II_WITH_P4EST
325namespace parallel
326{
327 namespace distributed
328 {
438 template <int dim, typename DataType>
440 {
441 public:
442 static_assert(
443 std::is_base_of_v<TransferableQuadraturePointData, DataType>,
444 "User's DataType class should be derived from TransferableQuadraturePointData");
445
451
463 const Quadrature<dim> &mass_quadrature,
464 const Quadrature<dim> &data_quadrature);
465
476 void
480
491 void
493
494 private:
500 std::vector<char>
503 &cell,
504 const CellStatus status);
505
511 void
514 &cell,
515 const CellStatus status,
516 const boost::iterator_range<std::vector<char>::const_iterator>
517 &data_range);
518
522 const std::unique_ptr<const FiniteElement<dim>> projection_fe;
523
529
533 const unsigned int n_q_points;
534
540
546
554
559
565
570 unsigned int handle;
571
576
582 };
583
584 } // namespace distributed
585
586} // namespace parallel
587
588#endif
589
592#ifndef DOXYGEN
593
594// ------------------- inline and template functions ----------------
595
596//--------------------------------------------------------------------
597// CellDataStorage
598//--------------------------------------------------------------------
599
600template <typename CellIteratorType, typename DataType>
601template <typename T>
602inline void
604 const CellIteratorType &cell,
605 const unsigned int n_q_points)
606{
607 static_assert(std::is_base_of_v<DataType, T>,
608 "User's T class should be derived from user's DataType class");
609 // The first time this method is called, it has to initialize the reference
610 // to the triangulation object
611 if (!tria)
612 tria = &cell->get_triangulation();
613 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
614
615 const auto key = cell->id();
616 if (map.find(key) == map.end())
617 {
618 map[key] = std::vector<std::shared_ptr<DataType>>(n_q_points);
619 // we need to initialize one-by-one as the std::vector<>(q, T())
620 // will end with a single same T object stored in each element of the
621 // vector:
622 const auto it = map.find(key);
623 for (unsigned int q = 0; q < n_q_points; ++q)
624 it->second[q] = std::make_shared<T>();
625 }
626}
627
628
629
630template <typename CellIteratorType, typename DataType>
631template <typename T>
632inline void
634 const CellIteratorType &cell_start,
636 const unsigned int number)
637{
638 for (CellIteratorType it = cell_start; it != cell_end; ++it)
639 if (it->is_locally_owned())
640 initialize<T>(it, number);
641}
642
643
644
645template <typename CellIteratorType, typename DataType>
646inline bool
647CellDataStorage<CellIteratorType, DataType>::erase(const CellIteratorType &cell)
648{
649 const auto key = cell->id();
650 const auto it = map.find(key);
651 if (it == map.end())
652 return false;
653 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
654 for (unsigned int i = 0; i < it->second.size(); ++i)
655 {
656 Assert(
657 it->second[i].unique(),
659 "Can not erase the cell data multiple objects reference its data."));
660 }
661
662 return (map.erase(key) == 1);
663}
664
665
666
667template <typename CellIteratorType, typename DataType>
668inline void
670{
671 // Do not call
672 // map.clear();
673 // as we want to be sure no one uses the stored objects. Loop manually:
674 auto it = map.begin();
675 while (it != map.end())
676 {
677 // loop over all objects and see if no one is using them
678 for (unsigned int i = 0; i < it->second.size(); ++i)
679 {
680 Assert(
681 it->second[i].unique(),
683 "Can not erase the cell data, multiple objects reference it."));
684 }
685 it = map.erase(it);
686 }
687}
688
689
690
691template <typename CellIteratorType, typename DataType>
692template <typename T>
693inline std::vector<std::shared_ptr<T>>
695 const CellIteratorType &cell)
696{
697 static_assert(std::is_base_of_v<DataType, T>,
698 "User's T class should be derived from user's DataType class");
699 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
700
701 const auto it = map.find(cell->id());
702 Assert(it != map.end(), ExcMessage("Could not find data for the cell"));
703
704 // It would be nice to have a specialized version of this function for
705 // T==DataType. However explicit (i.e full) specialization of a member
706 // template is only allowed when the enclosing class is also explicitly (i.e
707 // fully) specialized. Thus, stick with copying of shared pointers even when
708 // the T==DataType:
709 std::vector<std::shared_ptr<T>> res(it->second.size());
710 for (unsigned int q = 0; q < res.size(); ++q)
711 {
712 res[q] = std::dynamic_pointer_cast<T>(it->second[q]);
713 Assert(res[q], ExcCellDataTypeMismatch());
714 }
715 return res;
716}
717
718
719
720template <typename CellIteratorType, typename DataType>
721template <typename T>
722inline std::vector<std::shared_ptr<const T>>
724 const CellIteratorType &cell) const
725{
726 static_assert(std::is_base_of_v<DataType, T>,
727 "User's T class should be derived from user's DataType class");
728 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
729
730 const auto it = map.find(cell->id());
731 Assert(it != map.end(), ExcMessage("Could not find QP data for the cell"));
732
733 // Cast base class to the desired class. This has to be done irrespectively of
734 // T==DataType as we need to return shared_ptr<const T> to make sure the user
735 // does not modify the content of QP objects
736 std::vector<std::shared_ptr<const T>> res(it->second.size());
737 for (unsigned int q = 0; q < res.size(); ++q)
738 {
739 res[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
740 Assert(res[q], ExcCellDataTypeMismatch());
741 }
742 return res;
743}
744
745template <typename CellIteratorType, typename DataType>
746template <typename T>
747inline std::optional<std::vector<std::shared_ptr<T>>>
749 const CellIteratorType &cell)
750{
751 static_assert(std::is_base_of_v<DataType, T>,
752 "User's T class should be derived from user's DataType class");
753 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
754
755 const auto it = map.find(cell->id());
756 if (it != map.end())
757 {
758 // Cast base class to the desired class. This has to be done
759 // irrespectively of T==DataType as we need to return
760 // shared_ptr<const T> to make sure the user
761 // does not modify the content of QP objects
762 std::vector<std::shared_ptr<T>> result(it->second.size());
763 for (unsigned int q = 0; q < result.size(); ++q)
764 {
765 result[q] = std::dynamic_pointer_cast<T>(it->second[q]);
766 Assert(result[q], ExcCellDataTypeMismatch());
767 }
768 return {result};
769 }
770 else
771 {
772 return {};
773 }
774}
775
776template <typename CellIteratorType, typename DataType>
777template <typename T>
778inline std::optional<std::vector<std::shared_ptr<const T>>>
780 const CellIteratorType &cell) const
781{
782 static_assert(std::is_base_of_v<DataType, T>,
783 "User's T class should be derived from user's DataType class");
784 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
785
786 const auto it = map.find(cell->id());
787 if (it != map.end())
788 {
789 // Cast base class to the desired class. This has to be done
790 // irrespectively of T==DataType as we need to return
791 // shared_ptr<const T> to make sure the user
792 // does not modify the content of QP objects
793 std::vector<std::shared_ptr<const T>> result(it->second.size());
794 for (unsigned int q = 0; q < result.size(); ++q)
795 {
796 result[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
797 Assert(result[q], ExcCellDataTypeMismatch());
798 }
799 return {result};
800 }
801 else
802 {
803 return {};
804 }
805}
806
807//--------------------------------------------------------------------
808// ContinuousQuadratureDataTransfer
809//--------------------------------------------------------------------
810
811
812/*
813 * Pack cell data of type @p DataType stored using @p data_storage in @p cell
814 * at each quadrature point to @p matrix_data. Here @p matrix_data is a matrix
815 * whose first index corresponds to different quadrature points on the cell
816 * whereas the second index represents different values stored at each
817 * quadrature point in the DataType class.
818 */
819template <typename CellIteratorType, typename DataType>
820inline void
821pack_cell_data(const CellIteratorType &cell,
823 FullMatrix<double> &matrix_data)
824{
825 static_assert(std::is_base_of_v<TransferableQuadraturePointData, DataType>,
826 "User's DataType class should be derived from QPData");
827
828 if (const auto qpd = data_storage->try_get_data(cell))
829 {
830 const unsigned int m = qpd->size();
831 Assert(m > 0, ExcInternalError());
832 const unsigned int n = (*qpd)[0]->number_of_values();
833 matrix_data.reinit(m, n);
834
835 std::vector<double> single_qp_data(n);
836 for (unsigned int q = 0; q < m; ++q)
837 {
838 (*qpd)[q]->pack_values(single_qp_data);
839 AssertDimension(single_qp_data.size(), n);
840
841 for (unsigned int i = 0; i < n; ++i)
842 matrix_data(q, i) = single_qp_data[i];
843 }
844 }
845 else
846 {
847 matrix_data.reinit({0, 0});
848 }
849}
850
851
852
853/*
854 * the opposite of the pack function above.
855 */
856template <typename CellIteratorType, typename DataType>
857inline void
858unpack_to_cell_data(const CellIteratorType &cell,
859 const FullMatrix<double> &values_at_qp,
861{
862 static_assert(std::is_base_of_v<TransferableQuadraturePointData, DataType>,
863 "User's DataType class should be derived from QPData");
864
865 if (const auto qpd = data_storage->try_get_data(cell))
866 {
867 const unsigned int n = values_at_qp.n();
868 AssertDimension((*qpd)[0]->number_of_values(), n);
869
870 std::vector<double> single_qp_data(n);
871 AssertDimension(qpd->size(), values_at_qp.m());
872
873 for (unsigned int q = 0; q < qpd->size(); ++q)
874 {
875 for (unsigned int i = 0; i < n; ++i)
876 single_qp_data[i] = values_at_qp(q, i);
877 (*qpd)[q]->unpack_values(single_qp_data);
878 }
879 }
880}
881
882
883# ifdef DEAL_II_WITH_P4EST
884
885namespace parallel
886{
887 namespace distributed
888 {
889 template <int dim, typename DataType>
892 const Quadrature<dim> &lhs_quadrature,
893 const Quadrature<dim> &rhs_quadrature)
894 : projection_fe(
895 std::unique_ptr<const FiniteElement<dim>>(projection_fe_.clone()))
896 , data_size_in_bytes(0)
897 , n_q_points(rhs_quadrature.size())
898 , project_to_fe_matrix(projection_fe->n_dofs_per_cell(), n_q_points)
899 , project_to_qp_matrix(n_q_points, projection_fe->n_dofs_per_cell())
900 , handle(numbers::invalid_unsigned_int)
901 , data_storage(nullptr)
902 , triangulation(nullptr)
903 {
904 Assert(
905 projection_fe->n_components() == 1,
907 "ContinuousQuadratureDataTransfer requires scalar FiniteElement"));
908
910 *projection_fe.get(),
911 lhs_quadrature,
912 rhs_quadrature,
914
916 *projection_fe.get(), rhs_quadrature, project_to_qp_matrix);
917 }
918
919
920
921 template <int dim, typename DataType>
922 inline void
927 {
928 Assert(data_storage == nullptr,
929 ExcMessage("This function can be called only once"));
930 triangulation = &tr_;
931 data_storage = &data_storage_;
932
933 handle = triangulation->register_data_attach(
934 [this](const typename parallel::distributed::Triangulation<
935 dim>::cell_iterator &cell,
936 const CellStatus status) {
937 return this->pack_function(cell, status);
938 },
939 /*returns_variable_size_data=*/true);
940 }
941
942
943
944 template <int dim, typename DataType>
945 inline void
947 {
948 triangulation->notify_ready_to_unpack(
949 handle,
950 [this](const typename parallel::distributed::Triangulation<
951 dim>::cell_iterator &cell,
952 const CellStatus status,
953 const boost::iterator_range<std::vector<char>::const_iterator>
954 &data_range) {
955 this->unpack_function(cell, status, data_range);
956 });
957
958 // invalidate the pointers
959 data_storage = nullptr;
960 triangulation = nullptr;
961 }
962
963
964
965 template <int dim, typename DataType>
966 inline std::vector<char>
969 &cell,
970 const CellStatus /*status*/)
971 {
972 pack_cell_data(cell, data_storage, matrix_quadrature);
973
974 // project to FE
975 const unsigned int number_of_values = matrix_quadrature.n();
976 matrix_dofs.reinit(project_to_fe_matrix.m(), number_of_values);
977 if (number_of_values > 0)
978 project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
979
980 return Utilities::pack(matrix_dofs, /*allow_compression=*/false);
981 }
982
983
984
985 template <int dim, typename DataType>
986 inline void
989 &cell,
990 const CellStatus status,
991 const boost::iterator_range<std::vector<char>::const_iterator>
992 &data_range)
993 {
996 (void)status;
997
998 matrix_dofs =
999 Utilities::unpack<FullMatrix<double>>(data_range.begin(),
1000 data_range.end(),
1001 /*allow_compression=*/false);
1002 const unsigned int number_of_values = matrix_dofs.n();
1003 if (number_of_values == 0)
1004 return;
1005
1006 matrix_quadrature.reinit(n_q_points, number_of_values);
1007
1008 if (cell->has_children())
1009 {
1010 // we need to first use prolongation matrix to get dofvalues on child
1011 // cells based on dofvalues stored in the parent's data_store
1012 matrix_dofs_child.reinit(projection_fe->n_dofs_per_cell(),
1013 number_of_values);
1014 for (unsigned int child = 0; child < cell->n_children(); ++child)
1015 if (cell->child(child)->is_locally_owned())
1016 {
1017 projection_fe
1018 ->get_prolongation_matrix(child, cell->refinement_case())
1019 .mmult(matrix_dofs_child, matrix_dofs);
1020
1021 // now we do the usual business of evaluating FE on quadrature
1022 // points:
1023 project_to_qp_matrix.mmult(matrix_quadrature,
1024 matrix_dofs_child);
1025
1026 // finally, put back into the map:
1027 unpack_to_cell_data(cell->child(child),
1028 matrix_quadrature,
1029 data_storage);
1030 }
1031 }
1032 else
1033 {
1034 // if there are no children, evaluate FE field at
1035 // rhs_quadrature points.
1036 project_to_qp_matrix.mmult(matrix_quadrature, matrix_dofs);
1037
1038 // finally, put back into the map:
1039 unpack_to_cell_data(cell, matrix_quadrature, data_storage);
1040 }
1041 }
1042
1043 } // namespace distributed
1044
1045} // namespace parallel
1046
1047# endif // DEAL_II_WITH_P4EST
1048
1049#endif // DOXYGEN
1051
1052#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_start, const typename std_cxx20::type_identity< CellIteratorType >::type &cell_end, const unsigned int number_of_data_points_per_cell)
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
std::vector< std::shared_ptr< const T > > get_data(const CellIteratorType &cell) const
ObserverPointer< const Triangulation< dimension, space_dimension >, CellDataStorage< CellIteratorType, DataType > > tria
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
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcCellDataTypeMismatch()
#define AssertDimension(dim1, dim2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:489
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcTriangulationMismatch()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition tria.h:324
std::size_t size
Definition mpi.cc:734
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:1381
STL namespace.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation