Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14: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 - 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 
23 
25 
26 #include <deal.II/fe/fe.h>
27 #include <deal.II/fe/fe_tools.h>
28 
29 #include <deal.II/grid/tria.h>
31 
32 #include <deal.II/lac/vector.h>
33 
34 #include <map>
35 #include <optional>
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::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 
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_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 
528  std::size_t data_size_in_bytes;
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 
600 template <typename CellIteratorType, typename DataType>
601 template <typename T>
602 inline 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 
630 template <typename CellIteratorType, typename DataType>
631 template <typename T>
632 inline void
634  const CellIteratorType &cell_start,
635  const CellIteratorType &cell_end,
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 
645 template <typename CellIteratorType, typename DataType>
646 inline bool
647 CellDataStorage<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(),
658  ExcMessage(
659  "Can not erase the cell data multiple objects reference its data."));
660  }
661 
662  return (map.erase(key) == 1);
663 }
664 
665 
666 
667 template <typename CellIteratorType, typename DataType>
668 inline 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(),
682  ExcMessage(
683  "Can not erase the cell data, multiple objects reference it."));
684  }
685  it = map.erase(it);
686  }
687 }
688 
689 
690 
691 template <typename CellIteratorType, typename DataType>
692 template <typename T>
693 inline 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 
720 template <typename CellIteratorType, typename DataType>
721 template <typename T>
722 inline 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 
745 template <typename CellIteratorType, typename DataType>
746 template <typename T>
747 inline 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 
776 template <typename CellIteratorType, typename DataType>
777 template <typename T>
778 inline 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  */
819 template <typename CellIteratorType, typename DataType>
820 inline void
821 pack_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  */
856 template <typename CellIteratorType, typename DataType>
857 inline void
858 unpack_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 
885 namespace 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,
906  ExcMessage(
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:32
@ children_will_be_coarsened
std::vector< std::shared_ptr< const T > > get_data(const CellIteratorType &cell) const
CellDataStorage()=default
std::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
std::optional< std::vector< std::shared_ptr< const T > > > try_get_data(const CellIteratorType &cell) const
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)
~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
std::vector< char > pack_function(const typename parallel::distributed::Triangulation< dim >::cell_iterator &cell, const CellStatus status)
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 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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcTriangulationMismatch()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcCellDataTypeMismatch()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
static ::ExceptionBase & ExcTriangulationMismatch()
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:495
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:1343
static const unsigned int invalid_unsigned_int
Definition: types.h:221
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria