Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2016 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_quadrature_point_data_h
17#define dealii_quadrature_point_data_h
18
19#include <deal.II/base/config.h>
20
24
26
27#include <deal.II/fe/fe.h>
28#include <deal.II/fe/fe_tools.h>
29
30#include <deal.II/grid/tria.h>
32
33#include <deal.II/lac/vector.h>
34
35#include <map>
36#include <type_traits>
37#include <vector>
38
40
63template <typename CellIteratorType, typename DataType>
65{
66public:
70 CellDataStorage() = default;
71
75 ~CellDataStorage() override = default;
76
102 template <typename T = DataType>
103 void
104 initialize(const CellIteratorType &cell,
105 const unsigned int number_of_data_points_per_cell);
106
112 template <typename T = DataType>
113 void
114 initialize(const CellIteratorType &cell_start,
115 const CellIteratorType &cell_end,
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_cxx17::optional<std::vector<std::shared_ptr<T>>>
190 try_get_data(const CellIteratorType &cell);
191
208 template <typename T = DataType>
209 std_cxx17::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<TransferableQuadraturePointData, DataType>::value,
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,
505 status);
506
512 void
515 &cell,
517 status,
518 const boost::iterator_range<std::vector<char>::const_iterator>
519 &data_range);
520
524 const std::unique_ptr<const FiniteElement<dim>> projection_fe;
525
531
535 const unsigned int n_q_points;
536
542
548
556
561
567
572 unsigned int handle;
573
578
584 };
585
586 } // namespace distributed
587
588} // namespace parallel
589
590#endif
591
594#ifndef DOXYGEN
595
596// ------------------- inline and template functions ----------------
597
598//--------------------------------------------------------------------
599// CellDataStorage
600//--------------------------------------------------------------------
601
602template <typename CellIteratorType, typename DataType>
603template <typename T>
604inline void
606 const CellIteratorType &cell,
607 const unsigned int n_q_points)
608{
609 static_assert(std::is_base_of<DataType, T>::value,
610 "User's T class should be derived from user's DataType class");
611 // The first time this method is called, it has to initialize the reference
612 // to the triangulation object
613 if (!tria)
614 tria = &cell->get_triangulation();
615 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
616
617 const auto key = cell->id();
618 if (map.find(key) == map.end())
619 {
620 map[key] = std::vector<std::shared_ptr<DataType>>(n_q_points);
621 // we need to initialize one-by-one as the std::vector<>(q, T())
622 // will end with a single same T object stored in each element of the
623 // vector:
624 const auto it = map.find(key);
625 for (unsigned int q = 0; q < n_q_points; ++q)
626 it->second[q] = std::make_shared<T>();
627 }
628}
629
630
631
632template <typename CellIteratorType, typename DataType>
633template <typename T>
634inline void
636 const CellIteratorType &cell_start,
637 const CellIteratorType &cell_end,
638 const unsigned int number)
639{
640 for (CellIteratorType it = cell_start; it != cell_end; ++it)
641 if (it->is_locally_owned())
642 initialize<T>(it, number);
643}
644
645
646
647template <typename CellIteratorType, typename DataType>
648inline bool
649CellDataStorage<CellIteratorType, DataType>::erase(const CellIteratorType &cell)
650{
651 const auto key = cell->id();
652 const auto it = map.find(key);
653 if (it == map.end())
654 return false;
655 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
656 for (unsigned int i = 0; i < it->second.size(); ++i)
657 {
658 Assert(
659 it->second[i].unique(),
661 "Can not erase the cell data multiple objects reference its data."));
662 }
663
664 return (map.erase(key) == 1);
665}
666
667
668
669template <typename CellIteratorType, typename DataType>
670inline void
672{
673 // Do not call
674 // map.clear();
675 // as we want to be sure no one uses the stored objects. Loop manually:
676 auto it = map.begin();
677 while (it != map.end())
678 {
679 // loop over all objects and see if no one is using them
680 for (unsigned int i = 0; i < it->second.size(); ++i)
681 {
682 Assert(
683 it->second[i].unique(),
685 "Can not erase the cell data, multiple objects reference it."));
686 }
687 it = map.erase(it);
688 }
689}
690
691
692
693template <typename CellIteratorType, typename DataType>
694template <typename T>
695inline std::vector<std::shared_ptr<T>>
697 const CellIteratorType &cell)
698{
699 static_assert(std::is_base_of<DataType, T>::value,
700 "User's T class should be derived from user's DataType class");
701 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
702
703 const auto it = map.find(cell->id());
704 Assert(it != map.end(), ExcMessage("Could not find data for the cell"));
705
706 // It would be nice to have a specialized version of this function for
707 // T==DataType. However explicit (i.e full) specialization of a member
708 // template is only allowed when the enclosing class is also explicitly (i.e
709 // fully) specialized. Thus, stick with copying of shared pointers even when
710 // the T==DataType:
711 std::vector<std::shared_ptr<T>> res(it->second.size());
712 for (unsigned int q = 0; q < res.size(); ++q)
713 {
714 res[q] = std::dynamic_pointer_cast<T>(it->second[q]);
715 Assert(res[q], ExcCellDataTypeMismatch());
716 }
717 return res;
718}
719
720
721
722template <typename CellIteratorType, typename DataType>
723template <typename T>
724inline std::vector<std::shared_ptr<const T>>
726 const CellIteratorType &cell) const
727{
728 static_assert(std::is_base_of<DataType, T>::value,
729 "User's T class should be derived from user's DataType class");
730 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
731
732 const auto it = map.find(cell->id());
733 Assert(it != map.end(), ExcMessage("Could not find QP data for the cell"));
734
735 // Cast base class to the desired class. This has to be done irrespectively of
736 // T==DataType as we need to return shared_ptr<const T> to make sure the user
737 // does not modify the content of QP objects
738 std::vector<std::shared_ptr<const T>> res(it->second.size());
739 for (unsigned int q = 0; q < res.size(); ++q)
740 {
741 res[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
742 Assert(res[q], ExcCellDataTypeMismatch());
743 }
744 return res;
745}
746
747template <typename CellIteratorType, typename DataType>
748template <typename T>
749inline std_cxx17::optional<std::vector<std::shared_ptr<T>>>
751 const CellIteratorType &cell)
752{
753 static_assert(std::is_base_of<DataType, T>::value,
754 "User's T class should be derived from user's DataType class");
755 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
756
757 const auto it = map.find(cell->id());
758 if (it != map.end())
759 {
760 // Cast base class to the desired class. This has to be done
761 // irrespectively of T==DataType as we need to return
762 // shared_ptr<const T> to make sure the user
763 // does not modify the content of QP objects
764 std::vector<std::shared_ptr<T>> result(it->second.size());
765 for (unsigned int q = 0; q < result.size(); ++q)
766 {
767 result[q] = std::dynamic_pointer_cast<T>(it->second[q]);
768 Assert(result[q], ExcCellDataTypeMismatch());
769 }
770 return {result};
771 }
772 else
773 {
774 return {};
775 }
776}
777
778template <typename CellIteratorType, typename DataType>
779template <typename T>
780inline std_cxx17::optional<std::vector<std::shared_ptr<const T>>>
782 const CellIteratorType &cell) const
783{
784 static_assert(std::is_base_of<DataType, T>::value,
785 "User's T class should be derived from user's DataType class");
786 Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
787
788 const auto it = map.find(cell->id());
789 if (it != map.end())
790 {
791 // Cast base class to the desired class. This has to be done
792 // irrespectively of T==DataType as we need to return
793 // shared_ptr<const T> to make sure the user
794 // does not modify the content of QP objects
795 std::vector<std::shared_ptr<const T>> result(it->second.size());
796 for (unsigned int q = 0; q < result.size(); ++q)
797 {
798 result[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
799 Assert(result[q], ExcCellDataTypeMismatch());
800 }
801 return {result};
802 }
803 else
804 {
805 return {};
806 }
807}
808
809//--------------------------------------------------------------------
810// ContinuousQuadratureDataTransfer
811//--------------------------------------------------------------------
812
813
814/*
815 * Pack cell data of type @p DataType stored using @p data_storage in @p cell
816 * at each quadrature point to @p matrix_data. Here @p matrix_data is a matrix
817 * whose first index corresponds to different quadrature points on the cell
818 * whereas the second index represents different values stored at each
819 * quadrature point in the DataType class.
820 */
821template <typename CellIteratorType, typename DataType>
822inline void
823pack_cell_data(const CellIteratorType & cell,
825 FullMatrix<double> & matrix_data)
826{
827 static_assert(
828 std::is_base_of<TransferableQuadraturePointData, DataType>::value,
829 "User's DataType class should be derived from QPData");
830
831 if (const auto qpd = data_storage->try_get_data(cell))
832 {
833 const unsigned int m = qpd->size();
834 Assert(m > 0, ExcInternalError());
835 const unsigned int n = (*qpd)[0]->number_of_values();
836 matrix_data.reinit(m, n);
837
838 std::vector<double> single_qp_data(n);
839 for (unsigned int q = 0; q < m; ++q)
840 {
841 (*qpd)[q]->pack_values(single_qp_data);
842 AssertDimension(single_qp_data.size(), n);
843
844 for (unsigned int i = 0; i < n; ++i)
845 matrix_data(q, i) = single_qp_data[i];
846 }
847 }
848 else
849 {
850 matrix_data.reinit({0, 0});
851 }
852}
853
854
855
856/*
857 * the opposite of the pack function above.
858 */
859template <typename CellIteratorType, typename DataType>
860inline void
861unpack_to_cell_data(const CellIteratorType & cell,
862 const FullMatrix<double> & values_at_qp,
864{
865 static_assert(
866 std::is_base_of<TransferableQuadraturePointData, DataType>::value,
867 "User's DataType class should be derived from QPData");
868
869 if (const auto qpd = data_storage->try_get_data(cell))
870 {
871 const unsigned int n = values_at_qp.n();
872 AssertDimension((*qpd)[0]->number_of_values(), n);
873
874 std::vector<double> single_qp_data(n);
875 AssertDimension(qpd->size(), values_at_qp.m());
876
877 for (unsigned int q = 0; q < qpd->size(); ++q)
878 {
879 for (unsigned int i = 0; i < n; ++i)
880 single_qp_data[i] = values_at_qp(q, i);
881 (*qpd)[q]->unpack_values(single_qp_data);
882 }
883 }
884}
885
886
887# ifdef DEAL_II_WITH_P4EST
888
889namespace parallel
890{
891 namespace distributed
892 {
893 template <int dim, typename DataType>
896 const Quadrature<dim> & lhs_quadrature,
897 const Quadrature<dim> & rhs_quadrature)
898 : projection_fe(
899 std::unique_ptr<const FiniteElement<dim>>(projection_fe_.clone()))
900 , data_size_in_bytes(0)
901 , n_q_points(rhs_quadrature.size())
902 , project_to_fe_matrix(projection_fe->n_dofs_per_cell(), n_q_points)
903 , project_to_qp_matrix(n_q_points, projection_fe->n_dofs_per_cell())
904 , handle(numbers::invalid_unsigned_int)
905 , data_storage(nullptr)
906 , triangulation(nullptr)
907 {
908 Assert(
909 projection_fe->n_components() == 1,
911 "ContinuousQuadratureDataTransfer requires scalar FiniteElement"));
912
914 *projection_fe.get(),
915 lhs_quadrature,
916 rhs_quadrature,
918
920 *projection_fe.get(), rhs_quadrature, project_to_qp_matrix);
921 }
922
923
924
925 template <int dim, typename DataType>
926 inline void
931 {
932 Assert(data_storage == nullptr,
933 ExcMessage("This function can be called only once"));
934 triangulation = &tr_;
935 data_storage = &data_storage_;
936
937 handle = triangulation->register_data_attach(
938 [this](
940 dim>::cell_iterator &cell,
942 status) { return this->pack_function(cell, status); },
943 /*returns_variable_size_data=*/true);
944 }
945
946
947
948 template <int dim, typename DataType>
949 inline void
951 {
952 triangulation->notify_ready_to_unpack(
953 handle,
954 [this](
956 dim>::cell_iterator &cell,
958 status,
959 const boost::iterator_range<std::vector<char>::const_iterator>
960 &data_range) { this->unpack_function(cell, status, data_range); });
961
962 // invalidate the pointers
963 data_storage = nullptr;
964 triangulation = nullptr;
965 }
966
967
968
969 template <int dim, typename DataType>
970 inline std::vector<char>
973 &cell,
975 dim>::CellStatus /*status*/)
976 {
977 pack_cell_data(cell, data_storage, matrix_quadrature);
978
979 // project to FE
980 const unsigned int number_of_values = matrix_quadrature.n();
981 matrix_dofs.reinit(project_to_fe_matrix.m(), number_of_values);
982 if (number_of_values > 0)
983 project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
984
985 return Utilities::pack(matrix_dofs, /*allow_compression=*/false);
986 }
987
988
989
990 template <int dim, typename DataType>
991 inline void
994 &cell,
996 status,
997 const boost::iterator_range<std::vector<char>::const_iterator>
998 &data_range)
999 {
1000 Assert((status !=
1003 (void)status;
1004
1005 matrix_dofs =
1006 Utilities::unpack<FullMatrix<double>>(data_range.begin(),
1007 data_range.end(),
1008 /*allow_compression=*/false);
1009 const unsigned int number_of_values = matrix_dofs.n();
1010 if (number_of_values == 0)
1011 return;
1012
1013 matrix_quadrature.reinit(n_q_points, number_of_values);
1014
1015 if (cell->has_children())
1016 {
1017 // we need to first use prolongation matrix to get dofvalues on child
1018 // cells based on dofvalues stored in the parent's data_store
1019 matrix_dofs_child.reinit(projection_fe->n_dofs_per_cell(),
1020 number_of_values);
1021 for (unsigned int child = 0; child < cell->n_children(); ++child)
1022 if (cell->child(child)->is_locally_owned())
1023 {
1024 projection_fe
1025 ->get_prolongation_matrix(child, cell->refinement_case())
1026 .mmult(matrix_dofs_child, matrix_dofs);
1027
1028 // now we do the usual business of evaluating FE on quadrature
1029 // points:
1030 project_to_qp_matrix.mmult(matrix_quadrature,
1031 matrix_dofs_child);
1032
1033 // finally, put back into the map:
1034 unpack_to_cell_data(cell->child(child),
1035 matrix_quadrature,
1036 data_storage);
1037 }
1038 }
1039 else
1040 {
1041 // if there are no children, evaluate FE field at
1042 // rhs_quadrature points.
1043 project_to_qp_matrix.mmult(matrix_quadrature, matrix_dofs);
1044
1045 // finally, put back into the map:
1046 unpack_to_cell_data(cell, matrix_quadrature, data_storage);
1047 }
1048 }
1049
1050 } // namespace distributed
1051
1052} // namespace parallel
1053
1054# endif // DEAL_II_WITH_P4EST
1055
1056#endif // DOXYGEN
1058
1059#endif
std_cxx17::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
std_cxx17::optional< std::vector< std::shared_ptr< T > > > try_get_data(const CellIteratorType &cell)
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)
~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
void unpack_function(const typename parallel::distributed::Triangulation< dim >::cell_iterator &cell, const typename parallel::distributed::Triangulation< dim >::CellStatus status, const boost::iterator_range< std::vector< char >::const_iterator > &data_range)
void prepare_for_coarsening_and_refinement(parallel::distributed::Triangulation< dim > &tria, CellDataStorage< CellIteratorType, DataType > &data_storage)
std::vector< char > pack_function(const typename parallel::distributed::Triangulation< dim >::cell_iterator &cell, const typename parallel::distributed::Triangulation< dim >::CellStatus status)
CellDataStorage< CellIteratorType, DataType > * data_storage
typename ::Triangulation< dim, spacedim >::CellStatus CellStatus
Definition tria.h:294
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcCellDataTypeMismatch()
#define AssertDimension(dim1, dim2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcTriangulationMismatch()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition tria.h:270
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:1351
STL namespace.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria