Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 - 2019 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 
43 
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 
280 {
281 public:
286 
290  virtual ~TransferableQuadraturePointData() = default;
291 
297  virtual unsigned int
298  number_of_values() const = 0;
299 
310  virtual void
311  pack_values(std::vector<double> &values) const = 0;
312 
321  virtual void
322  unpack_values(const std::vector<double> &values) = 0;
323 };
324 
325 
326 #ifdef DEAL_II_WITH_P4EST
327 namespace parallel
328 {
329  namespace distributed
330  {
442  template <int dim, typename DataType>
444  {
445  public:
446  static_assert(
448  "User's DataType class should be derived from TransferableQuadraturePointData");
449 
453  using CellIteratorType =
455 
467  const Quadrature<dim> &mass_quadrature,
468  const Quadrature<dim> &data_quadrature);
469 
480  void
481  prepare_for_coarsening_and_refinement(
484 
495  void
496  interpolate();
497 
498  private:
504  std::vector<char>
505  pack_function(
507  &cell,
509  status);
510 
516  void
517  unpack_function(
519  &cell,
521  status,
522  const boost::iterator_range<std::vector<char>::const_iterator>
523  &data_range);
524 
528  const std::unique_ptr<const FiniteElement<dim>> projection_fe;
529 
534  std::size_t data_size_in_bytes;
535 
539  const unsigned int n_q_points;
540 
546 
552 
560 
565 
571 
576  unsigned int handle;
577 
582 
588  };
589 
590  } // namespace distributed
591 
592 } // namespace parallel
593 
594 #endif
595 
598 #ifndef DOXYGEN
599 
600 // ------------------- inline and template functions ----------------
601 
602 //--------------------------------------------------------------------
603 // CellDataStorage
604 //--------------------------------------------------------------------
605 
606 template <typename CellIteratorType, typename DataType>
607 template <typename T>
608 inline void
610  const CellIteratorType &cell,
611  const unsigned int n_q_points)
612 {
614  "User's T class should be derived from user's DataType class");
615  // The first time this method is called, it has to initialize the reference
616  // to the triangulation object
617  if (!tria)
618  tria = &cell->get_triangulation();
619  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
620 
621  const auto key = cell->id();
622  if (map.find(key) == map.end())
623  {
624  map[key] = std::vector<std::shared_ptr<DataType>>(n_q_points);
625  // we need to initialize one-by-one as the std::vector<>(q, T())
626  // will end with a single same T object stored in each element of the
627  // vector:
628  const auto it = map.find(key);
629  for (unsigned int q = 0; q < n_q_points; q++)
630  it->second[q] = std::make_shared<T>();
631  }
632 }
633 
634 
635 
636 template <typename CellIteratorType, typename DataType>
637 template <typename T>
638 inline void
640  const CellIteratorType &cell_start,
641  const CellIteratorType &cell_end,
642  const unsigned int number)
643 {
644  for (CellIteratorType it = cell_start; it != cell_end; it++)
645  if (it->is_locally_owned())
646  initialize<T>(it, number);
647 }
648 
649 
650 
651 template <typename CellIteratorType, typename DataType>
652 inline bool
654 {
655  const auto key = cell->id();
656  const auto it = map.find(key);
657  if (it == map.end())
658  return false;
659  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
660  for (unsigned int i = 0; i < it->second.size(); i++)
661  {
662  Assert(
663  it->second[i].unique(),
664  ExcMessage(
665  "Can not erase the cell data multiple objects reference its data."));
666  }
667 
668  return (map.erase(key) == 1);
669 }
670 
671 
672 
673 template <typename CellIteratorType, typename DataType>
674 inline void
676 {
677  // Do not call
678  // map.clear();
679  // as we want to be sure no one uses the stored objects. Loop manually:
680  auto it = map.begin();
681  while (it != map.end())
682  {
683  // loop over all objects and see if no one is using them
684  for (unsigned int i = 0; i < it->second.size(); i++)
685  {
686  Assert(
687  it->second[i].unique(),
688  ExcMessage(
689  "Can not erase the cell data, multiple objects reference it."));
690  }
691  it = map.erase(it);
692  }
693 }
694 
695 
696 
697 template <typename CellIteratorType, typename DataType>
698 template <typename T>
699 inline std::vector<std::shared_ptr<T>>
701  const CellIteratorType &cell)
702 {
704  "User's T class should be derived from user's DataType class");
705  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
706 
707  const auto it = map.find(cell->id());
708  Assert(it != map.end(), ExcMessage("Could not find data for the cell"));
709 
710  // It would be nice to have a specialized version of this function for
711  // T==DataType. However explicit (i.e full) specialization of a member
712  // template is only allowed when the enclosing class is also explicitly (i.e
713  // fully) specialized. Thus, stick with copying of shared pointers even when
714  // the T==DataType:
715  std::vector<std::shared_ptr<T>> res(it->second.size());
716  for (unsigned int q = 0; q < res.size(); q++)
717  {
718  res[q] = std::dynamic_pointer_cast<T>(it->second[q]);
719  Assert(res[q], ExcCellDataTypeMismatch());
720  }
721  return res;
722 }
723 
724 
725 
726 template <typename CellIteratorType, typename DataType>
727 template <typename T>
728 inline std::vector<std::shared_ptr<const T>>
730  const CellIteratorType &cell) const
731 {
733  "User's T class should be derived from user's DataType class");
734  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
735 
736  const auto it = map.find(cell->id());
737  Assert(it != map.end(), ExcMessage("Could not find QP data for the cell"));
738 
739  // Cast base class to the desired class. This has to be done irrespectively of
740  // T==DataType as we need to return shared_ptr<const T> to make sure the user
741  // does not modify the content of QP objects
742  std::vector<std::shared_ptr<const T>> res(it->second.size());
743  for (unsigned int q = 0; q < res.size(); q++)
744  {
745  res[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
746  Assert(res[q], ExcCellDataTypeMismatch());
747  }
748  return res;
749 }
750 
751 template <typename CellIteratorType, typename DataType>
752 template <typename T>
753 inline std_cxx17::optional<std::vector<std::shared_ptr<T>>>
755  const CellIteratorType &cell)
756 {
758  "User's T class should be derived from user's DataType class");
759  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
760 
761  const auto it = map.find(cell->id());
762  if (it != map.end())
763  {
764  // Cast base class to the desired class. This has to be done
765  // irrespectively of T==DataType as we need to return
766  // shared_ptr<const T> to make sure the user
767  // does not modify the content of QP objects
768  std::vector<std::shared_ptr<T>> result(it->second.size());
769  for (unsigned int q = 0; q < result.size(); q++)
770  {
771  result[q] = std::dynamic_pointer_cast<T>(it->second[q]);
772  Assert(result[q], ExcCellDataTypeMismatch());
773  }
774  return {result};
775  }
776  else
777  {
778  return {};
779  }
780 }
781 
782 template <typename CellIteratorType, typename DataType>
783 template <typename T>
784 inline std_cxx17::optional<std::vector<std::shared_ptr<const T>>>
786  const CellIteratorType &cell) const
787 {
789  "User's T class should be derived from user's DataType class");
790  Assert(&cell->get_triangulation() == tria, ExcTriangulationMismatch());
791 
792  const auto it = map.find(cell->id());
793  if (it != map.end())
794  {
795  // Cast base class to the desired class. This has to be done
796  // irrespectively of T==DataType as we need to return
797  // shared_ptr<const T> to make sure the user
798  // does not modify the content of QP objects
799  std::vector<std::shared_ptr<const T>> result(it->second.size());
800  for (unsigned int q = 0; q < result.size(); q++)
801  {
802  result[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
803  Assert(result[q], ExcCellDataTypeMismatch());
804  }
805  return {result};
806  }
807  else
808  {
809  return {};
810  }
811 }
812 
813 //--------------------------------------------------------------------
814 // ContinuousQuadratureDataTransfer
815 //--------------------------------------------------------------------
816 
817 
818 /*
819  * Pack cell data of type @p DataType stored using @p data_storage in @p cell
820  * at each quadrature point to @p matrix_data. Here @p matrix_data is a matrix
821  * whose first index corresponds to different quadrature points on the cell
822  * whereas the second index represents different values stored at each
823  * quadrature point in the DataType class.
824  */
825 template <typename CellIteratorType, typename DataType>
826 inline void
827 pack_cell_data(const CellIteratorType & cell,
829  FullMatrix<double> & matrix_data)
830 {
831  static_assert(
833  "User's DataType class should be derived from QPData");
834 
835  if (const auto qpd = data_storage->try_get_data(cell))
836  {
837  const unsigned int m = qpd->size();
838  Assert(m > 0, ExcInternalError());
839  const unsigned int n = (*qpd)[0]->number_of_values();
840  matrix_data.reinit(m, n);
841 
842  std::vector<double> single_qp_data(n);
843  for (unsigned int q = 0; q < m; ++q)
844  {
845  (*qpd)[q]->pack_values(single_qp_data);
846  AssertDimension(single_qp_data.size(), n);
847 
848  for (unsigned int i = 0; i < n; ++i)
849  matrix_data(q, i) = single_qp_data[i];
850  }
851  }
852  else
853  {
854  matrix_data.reinit({0, 0});
855  }
856 }
857 
858 
859 
860 /*
861  * the opposite of the pack function above.
862  */
863 template <typename CellIteratorType, typename DataType>
864 inline void
865 unpack_to_cell_data(const CellIteratorType & cell,
866  const FullMatrix<double> & values_at_qp,
868 {
869  static_assert(
871  "User's DataType class should be derived from QPData");
872 
873  if (const auto qpd = data_storage->try_get_data(cell))
874  {
875  const unsigned int n = values_at_qp.n();
876  AssertDimension((*qpd)[0]->number_of_values(), n);
877 
878  std::vector<double> single_qp_data(n);
879  AssertDimension(qpd->size(), values_at_qp.m());
880 
881  for (unsigned int q = 0; q < qpd->size(); ++q)
882  {
883  for (unsigned int i = 0; i < n; ++i)
884  single_qp_data[i] = values_at_qp(q, i);
885  (*qpd)[q]->unpack_values(single_qp_data);
886  }
887  }
888 }
889 
890 
891 # ifdef DEAL_II_WITH_P4EST
892 
893 namespace parallel
894 {
895  namespace distributed
896  {
897  template <int dim, typename DataType>
900  const Quadrature<dim> & lhs_quadrature,
901  const Quadrature<dim> & rhs_quadrature)
902  : projection_fe(
903  std::unique_ptr<const FiniteElement<dim>>(projection_fe_.clone()))
904  , data_size_in_bytes(0)
905  , n_q_points(rhs_quadrature.size())
906  , project_to_fe_matrix(projection_fe->dofs_per_cell, n_q_points)
907  , project_to_qp_matrix(n_q_points, projection_fe->dofs_per_cell)
909  , data_storage(nullptr)
910  , triangulation(nullptr)
911  {
912  Assert(
913  projection_fe->n_components() == 1,
914  ExcMessage(
915  "ContinuousQuadratureDataTransfer requires scalar FiniteElement"));
916 
918  *projection_fe.get(),
919  lhs_quadrature,
920  rhs_quadrature,
921  project_to_fe_matrix);
922 
924  *projection_fe.get(), rhs_quadrature, project_to_qp_matrix);
925  }
926 
927 
928 
929  template <int dim, typename DataType>
930  inline void
935  {
936  Assert(data_storage == nullptr,
937  ExcMessage("This function can be called only once"));
938  triangulation = &tr_;
939  data_storage = &data_storage_;
940 
941  handle = triangulation->register_data_attach(
942  [this](
944  dim>::cell_iterator &cell,
946  status) { return this->pack_function(cell, status); },
947  /*returns_variable_size_data=*/true);
948  }
949 
950 
951 
952  template <int dim, typename DataType>
953  inline void
955  {
956  triangulation->notify_ready_to_unpack(
957  handle,
958  [this](
960  dim>::cell_iterator &cell,
962  status,
963  const boost::iterator_range<std::vector<char>::const_iterator>
964  &data_range) { this->unpack_function(cell, status, data_range); });
965 
966  // invalidate the pointers
967  data_storage = nullptr;
968  triangulation = nullptr;
969  }
970 
971 
972 
973  template <int dim, typename DataType>
974  inline std::vector<char>
977  &cell,
979  dim>::CellStatus /*status*/)
980  {
981  pack_cell_data(cell, data_storage, matrix_quadrature);
982 
983  // project to FE
984  const unsigned int number_of_values = matrix_quadrature.n();
985  matrix_dofs.reinit(project_to_fe_matrix.m(), number_of_values);
986  if (number_of_values > 0)
987  project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
988 
989  return Utilities::pack(matrix_dofs, /*allow_compression=*/false);
990  }
991 
992 
993 
994  template <int dim, typename DataType>
995  inline void
998  &cell,
1000  status,
1001  const boost::iterator_range<std::vector<char>::const_iterator>
1002  &data_range)
1003  {
1004  Assert((status !=
1006  ExcNotImplemented());
1007  (void)status;
1008 
1009  matrix_dofs =
1010  Utilities::unpack<FullMatrix<double>>(data_range.begin(),
1011  data_range.end(),
1012  /*allow_compression=*/false);
1013  const unsigned int number_of_values = matrix_dofs.n();
1014  if (number_of_values == 0)
1015  return;
1016 
1017  matrix_quadrature.reinit(n_q_points, number_of_values);
1018 
1019  if (cell->has_children())
1020  {
1021  // we need to first use prolongation matrix to get dofvalues on child
1022  // cells based on dofvalues stored in the parent's data_store
1023  matrix_dofs_child.reinit(projection_fe->dofs_per_cell,
1024  number_of_values);
1025  for (unsigned int child = 0; child < cell->n_children(); ++child)
1026  if (cell->child(child)->is_locally_owned())
1027  {
1028  projection_fe
1029  ->get_prolongation_matrix(child, cell->refinement_case())
1030  .mmult(matrix_dofs_child, matrix_dofs);
1031 
1032  // now we do the usual business of evaluating FE on quadrature
1033  // points:
1034  project_to_qp_matrix.mmult(matrix_quadrature,
1035  matrix_dofs_child);
1036 
1037  // finally, put back into the map:
1038  unpack_to_cell_data(cell->child(child),
1039  matrix_quadrature,
1040  data_storage);
1041  }
1042  }
1043  else
1044  {
1045  // if there are no children, evaluate FE field at
1046  // rhs_quadrature points.
1047  project_to_qp_matrix.mmult(matrix_quadrature, matrix_dofs);
1048 
1049  // finally, put back into the map:
1050  unpack_to_cell_data(cell, matrix_quadrature, data_storage);
1051  }
1052  }
1053 
1054  } // namespace distributed
1055 
1056 } // namespace parallel
1057 
1058 # endif // DEAL_II_WITH_P4EST
1059 
1060 #endif // DOXYGEN
1062 
1063 #endif
SmartPointer< const Triangulation< dimension, space_dimension >, CellDataStorage< CellIteratorType, DataType > > tria
size_type m() const
~CellDataStorage() override=default
static const unsigned int invalid_unsigned_int
Definition: types.h:191
void compute_interpolation_to_quadrature_points_matrix(const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, FullMatrix< double > &I_q)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
CellDataStorage()=default
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const =0
std::vector< std::shared_ptr< T > > get_data(const CellIteratorType &cell)
bool erase(const CellIteratorType &cell)
CellDataStorage< CellIteratorType, DataType > * data_storage
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
static const char T
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1419
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
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)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
const std::unique_ptr< const FiniteElement< dim > > projection_fe
std::map< CellId, std::vector< std::shared_ptr< DataType > > > map
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1209
static ::ExceptionBase & ExcTriangulationMismatch()
unsigned int size() const
std_cxx17::optional< std::vector< std::shared_ptr< T > > > try_get_data(const CellIteratorType &cell)
static constexpr unsigned int dimension
static ::ExceptionBase & ExcCellDataTypeMismatch()
parallel::distributed::Triangulation< dim > * triangulation
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:361
static constexpr unsigned int space_dimension
static ::ExceptionBase & ExcNotImplemented()
typename parallel::distributed::Triangulation< dim >::cell_iterator CellIteratorType
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
static const bool value
void initialize(const CellIteratorType &cell, const unsigned int number_of_data_points_per_cell)
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
static ::ExceptionBase & ExcInternalError()
Triangulation< dim, spacedim > & get_triangulation()
Definition: tria.cc:13293