Reference documentation for deal.II version GIT a3f17f8a20 2023-06-02 10:50: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\}}\)
quadrature_point_data.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2021 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 
63 template <typename CellIteratorType, typename DataType>
65 {
66 public:
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
134  clear();
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 
212 private:
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 {
279 public:
284 
288  virtual ~TransferableQuadraturePointData() = default;
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
325 namespace 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 
530  std::size_t data_size_in_bytes;
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 
602 template <typename CellIteratorType, typename DataType>
603 template <typename T>
604 inline 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 
632 template <typename CellIteratorType, typename DataType>
633 template <typename T>
634 inline 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 
647 template <typename CellIteratorType, typename DataType>
648 inline bool
649 CellDataStorage<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(),
660  ExcMessage(
661  "Can not erase the cell data multiple objects reference its data."));
662  }
663 
664  return (map.erase(key) == 1);
665 }
666 
667 
668 
669 template <typename CellIteratorType, typename DataType>
670 inline 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(),
684  ExcMessage(
685  "Can not erase the cell data, multiple objects reference it."));
686  }
687  it = map.erase(it);
688  }
689 }
690 
691 
692 
693 template <typename CellIteratorType, typename DataType>
694 template <typename T>
695 inline 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 
722 template <typename CellIteratorType, typename DataType>
723 template <typename T>
724 inline 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 
747 template <typename CellIteratorType, typename DataType>
748 template <typename T>
749 inline 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 
778 template <typename CellIteratorType, typename DataType>
779 template <typename T>
780 inline 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  */
821 template <typename CellIteratorType, typename DataType>
822 inline void
823 pack_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  */
859 template <typename CellIteratorType, typename DataType>
860 inline void
861 unpack_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 
889 namespace 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,
910  ExcMessage(
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 !=
1002  ExcNotImplemented());
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::vector< std::shared_ptr< const T > > get_data(const CellIteratorType &cell) const
CellDataStorage()=default
std_cxx17::optional< std::vector< std::shared_ptr< T > > > try_get_data(const CellIteratorType &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
SmartPointer< const Triangulation< dimension, space_dimension >, CellDataStorage< CellIteratorType, DataType > > tria
std::vector< std::shared_ptr< T > > get_data(const CellIteratorType &cell)
void initialize(const CellIteratorType &cell_start, const CellIteratorType &cell_end, const unsigned int number_of_data_points_per_cell)
std_cxx17::optional< std::vector< std::shared_ptr< const T > > > try_get_data(const CellIteratorType &cell) const
~CellDataStorage() override=default
std::map< CellId, std::vector< std::shared_ptr< DataType > > > map
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)
std::vector< char > pack_function(const typename parallel::distributed::Triangulation< dim >::cell_iterator &cell, const typename parallel::distributed::Triangulation< dim >::CellStatus status)
void prepare_for_coarsening_and_refinement(parallel::distributed::Triangulation< dim > &tria, CellDataStorage< CellIteratorType, DataType > &data_storage)
CellDataStorage< CellIteratorType, DataType > * data_storage
typename ::Triangulation< dim, spacedim >::CellStatus CellStatus
Definition: tria.h:294
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcTriangulationMismatch()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcCellDataTypeMismatch()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
static ::ExceptionBase & ExcTriangulationMismatch()
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:488
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
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria