Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
quadrature_point_data.h
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 
21 #include <deal.II/base/quadrature.h>
22 #include <deal.II/base/subscriptor.h>
23 
24 #include <deal.II/distributed/tria.h>
25 
26 #include <deal.II/fe/fe.h>
27 #include <deal.II/fe/fe_tools.h>
28 
29 #include <deal.II/grid/tria.h>
30 #include <deal.II/grid/tria_accessor.h>
31 
32 #include <deal.II/lac/vector.h>
33 
34 #include <map>
35 #include <type_traits>
36 #include <vector>
37 
38 DEAL_II_NAMESPACE_OPEN
39 
42 
58 template <typename CellIteratorType, typename DataType>
60 {
61 public:
65  CellDataStorage() = default;
66 
70  ~CellDataStorage() override = default;
71 
93  template <typename T = DataType>
94  void
95  initialize(const CellIteratorType &cell,
96  const unsigned int number_of_data_points_per_cell);
97 
103  template <typename T = DataType>
104  void
105  initialize(const CellIteratorType &cell_start,
106  const CellIteratorType &cell_end,
107  const unsigned int number_of_data_points_per_cell);
108 
118  bool
119  erase(const CellIteratorType &cell);
120 
124  void
125  clear();
126 
138  template <typename T = DataType>
139  std::vector<std::shared_ptr<T>>
140  get_data(const CellIteratorType &cell);
141 
153  template <typename T = DataType>
154  std::vector<std::shared_ptr<const T>>
155  get_data(const CellIteratorType &cell) const;
156 
157 private:
161  std::map<CellIteratorType, std::vector<std::shared_ptr<DataType>>> map;
162 };
163 
164 
184 {
185 public:
190 
194  virtual ~TransferableQuadraturePointData() = default;
195 
201  virtual unsigned int
202  number_of_values() const = 0;
203 
214  virtual void
215  pack_values(std::vector<double> &values) const = 0;
216 
225  virtual void
226  unpack_values(const std::vector<double> &values) = 0;
227 };
228 
229 
230 #ifdef DEAL_II_WITH_P4EST
231 namespace parallel
232 {
233  namespace distributed
234  {
346  template <int dim, typename DataType>
348  {
349  public:
350  static_assert(
351  std::is_base_of<TransferableQuadraturePointData, DataType>::value,
352  "User's DataType class should be derived from TransferableQuadraturePointData");
353 
357  using CellIteratorType =
359 
371  const Quadrature<dim> &mass_quadrature,
372  const Quadrature<dim> &data_quadrature);
373 
384  void
388 
399  void
400  interpolate();
401 
402  private:
408  std::vector<char>
411  &cell,
412  const typename parallel::distributed::Triangulation<dim>::CellStatus
413  status);
414 
420  void
423  &cell,
424  const typename parallel::distributed::Triangulation<dim>::CellStatus
425  status,
426  const boost::iterator_range<std::vector<char>::const_iterator>
427  &data_range);
428 
432  const std::unique_ptr<const FiniteElement<dim>> projection_fe;
433 
438  std::size_t data_size_in_bytes;
439 
443  const unsigned int n_q_points;
444 
450 
456 
464 
469 
475 
480  unsigned int handle;
481 
486 
492  };
493 
494  } // namespace distributed
495 
496 } // namespace parallel
497 
498 #endif
499 
502 #ifndef DOXYGEN
503 
504 // ------------------- inline and template functions ----------------
505 
506 //--------------------------------------------------------------------
507 // CellDataStorage
508 //--------------------------------------------------------------------
509 
510 template <typename CellIteratorType, typename DataType>
511 template <typename T>
512 void
514  const CellIteratorType &cell,
515  const unsigned int n_q_points)
516 {
517  static_assert(std::is_base_of<DataType, T>::value,
518  "User's T class should be derived from user's DataType class");
519 
520  if (map.find(cell) == map.end())
521  {
522  map[cell] = std::vector<std::shared_ptr<DataType>>(n_q_points);
523  // we need to initialize one-by-one as the std::vector<>(q, T())
524  // will end with a single same T object stored in each element of the
525  // vector:
526  auto it = map.find(cell);
527  for (unsigned int q = 0; q < n_q_points; q++)
528  it->second[q] = std::make_shared<T>();
529  }
530 }
531 
532 
533 
534 template <typename CellIteratorType, typename DataType>
535 template <typename T>
536 void
538  const CellIteratorType &cell_start,
539  const CellIteratorType &cell_end,
540  const unsigned int number)
541 {
542  for (CellIteratorType it = cell_start; it != cell_end; it++)
543  if (it->is_locally_owned())
544  initialize<T>(it, number);
545 }
546 
547 
548 
549 template <typename CellIteratorType, typename DataType>
550 bool
551 CellDataStorage<CellIteratorType, DataType>::erase(const CellIteratorType &cell)
552 {
553  const auto it = map.find(cell);
554  for (unsigned int i = 0; i < it->second.size(); i++)
555  {
556  Assert(
557  it->second[i].unique(),
558  ExcMessage(
559  "Can not erase the cell data multiple objects reference its data."));
560  }
561 
562  return (map.erase(cell) == 1);
563 }
564 
565 
566 
567 template <typename CellIteratorType, typename DataType>
568 void
570 {
571  // Do not call
572  // map.clear();
573  // as we want to be sure no one uses the stored objects. Loop manually:
574  auto it = map.begin();
575  while (it != map.end())
576  {
577  // loop over all objects and see if no one is using them
578  for (unsigned int i = 0; i < it->second.size(); i++)
579  {
580  Assert(
581  it->second[i].unique(),
582  ExcMessage(
583  "Can not erase the cell data, multiple objects reference it."));
584  }
585  it = map.erase(it);
586  }
587 }
588 
589 
590 
591 template <typename CellIteratorType, typename DataType>
592 template <typename T>
593 std::vector<std::shared_ptr<T>>
595  const CellIteratorType &cell)
596 {
597  static_assert(std::is_base_of<DataType, T>::value,
598  "User's T class should be derived from user's DataType class");
599 
600  auto it = map.find(cell);
601  Assert(it != map.end(), ExcMessage("Could not find data for the cell"));
602 
603  // It would be nice to have a specialized version of this function for
604  // T==DataType. However explicit (i.e full) specialization of a member
605  // template is only allowed when the enclosing class is also explicitly (i.e
606  // fully) specialized. Thus, stick with copying of shared pointers even when
607  // the T==DataType:
608  std::vector<std::shared_ptr<T>> res(it->second.size());
609  for (unsigned int q = 0; q < res.size(); q++)
610  res[q] = std::dynamic_pointer_cast<T>(it->second[q]);
611  return res;
612 }
613 
614 
615 
616 template <typename CellIteratorType, typename DataType>
617 template <typename T>
618 std::vector<std::shared_ptr<const T>>
620  const CellIteratorType &cell) const
621 {
622  static_assert(std::is_base_of<DataType, T>::value,
623  "User's T class should be derived from user's DataType class");
624 
625  auto it = map.find(cell);
626  Assert(it != map.end(), ExcMessage("Could not find QP data for the cell"));
627 
628  // Cast base class to the desired class. This has to be done irrespectively of
629  // T==DataType as we need to return shared_ptr<const T> to make sure the user
630  // does not modify the content of QP objects
631  std::vector<std::shared_ptr<const T>> res(it->second.size());
632  for (unsigned int q = 0; q < res.size(); q++)
633  res[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
634  return res;
635 }
636 
637 //--------------------------------------------------------------------
638 // ContinuousQuadratureDataTransfer
639 //--------------------------------------------------------------------
640 
641 
642 /*
643  * Pack cell data of type @p DataType stored using @p data_storage in @p cell
644  * at each quadrature point to @p matrix_data. Here @p matrix_data is a matrix
645  * whose first index corresponds to different quadrature points on the cell
646  * whereas the second index represents different values stored at each
647  * quadrature point in the DataType class.
648  */
649 template <typename CellIteratorType, typename DataType>
650 void
651 pack_cell_data(const CellIteratorType & cell,
653  FullMatrix<double> & matrix_data)
654 {
655  static_assert(
656  std::is_base_of<TransferableQuadraturePointData, DataType>::value,
657  "User's DataType class should be derived from QPData");
658 
659  const std::vector<std::shared_ptr<const DataType>> qpd =
660  data_storage->get_data(cell);
661 
662  const unsigned int n = matrix_data.n();
663 
664  std::vector<double> single_qp_data(n);
665  Assert(qpd.size() == matrix_data.m(),
666  ExcDimensionMismatch(qpd.size(), matrix_data.m()));
667  for (unsigned int q = 0; q < qpd.size(); q++)
668  {
669  qpd[q]->pack_values(single_qp_data);
670  Assert(single_qp_data.size() == n,
671  ExcDimensionMismatch(single_qp_data.size(), n));
672 
673  for (unsigned int i = 0; i < n; i++)
674  matrix_data(q, i) = single_qp_data[i];
675  }
676 }
677 
678 
679 
680 /*
681  * the opposite of the pack function above.
682  */
683 template <typename CellIteratorType, typename DataType>
684 void
685 unpack_to_cell_data(const CellIteratorType & cell,
686  const FullMatrix<double> & values_at_qp,
688 {
689  static_assert(
690  std::is_base_of<TransferableQuadraturePointData, DataType>::value,
691  "User's DataType class should be derived from QPData");
692 
693  std::vector<std::shared_ptr<DataType>> qpd = data_storage->get_data(cell);
694 
695  const unsigned int n = values_at_qp.n();
696 
697  std::vector<double> single_qp_data(n);
698  Assert(qpd.size() == values_at_qp.m(),
699  ExcDimensionMismatch(qpd.size(), values_at_qp.m()));
700 
701  for (unsigned int q = 0; q < qpd.size(); q++)
702  {
703  for (unsigned int i = 0; i < n; i++)
704  single_qp_data[i] = values_at_qp(q, i);
705  qpd[q]->unpack_values(single_qp_data);
706  }
707 }
708 
709 
710 # ifdef DEAL_II_WITH_P4EST
711 
712 namespace parallel
713 {
714  namespace distributed
715  {
716  template <int dim, typename DataType>
719  const Quadrature<dim> & lhs_quadrature,
720  const Quadrature<dim> & rhs_quadrature)
721  : projection_fe(
722  std::unique_ptr<const FiniteElement<dim>>(projection_fe_.clone()))
723  , data_size_in_bytes(0)
724  , n_q_points(rhs_quadrature.size())
725  , project_to_fe_matrix(projection_fe->dofs_per_cell, n_q_points)
726  , project_to_qp_matrix(n_q_points, projection_fe->dofs_per_cell)
727  , handle(numbers::invalid_unsigned_int)
728  , data_storage(nullptr)
729  , triangulation(nullptr)
730  {
731  Assert(
732  projection_fe->n_components() == 1,
733  ExcMessage(
734  "ContinuousQuadratureDataTransfer requires scalar FiniteElement"));
735 
737  *projection_fe.get(),
738  lhs_quadrature,
739  rhs_quadrature,
741 
743  *projection_fe.get(), rhs_quadrature, project_to_qp_matrix);
744  }
745 
746 
747 
748  template <int dim, typename DataType>
749  void
754  {
755  Assert(data_storage == nullptr,
756  ExcMessage("This function can be called only once"));
757  triangulation = &tr_;
758  data_storage = &data_storage_;
759  // get the number from the first active cell
760  unsigned int number_of_values = 0;
761  // if triangulation has some active cells locally owned cells on this
762  // processor we can expect data to be initialized. Do that to get the
763  // number:
765  dim>::active_cell_iterator it = triangulation->begin_active();
766  it != triangulation->end();
767  it++)
768  if (it->is_locally_owned())
769  {
770  std::vector<std::shared_ptr<DataType>> qpd =
771  data_storage->get_data(it);
772  number_of_values = qpd[0]->number_of_values();
773  break;
774  }
775  // some processors may have no data stored, thus get the maximum among all
776  // processors:
777  number_of_values = Utilities::MPI::max(number_of_values,
778  triangulation->get_communicator());
779  Assert(number_of_values > 0, ExcInternalError());
780  const unsigned int dofs_per_cell = projection_fe->dofs_per_cell;
781  matrix_dofs.reinit(dofs_per_cell, number_of_values);
782  matrix_dofs_child.reinit(dofs_per_cell, number_of_values);
783  matrix_quadrature.reinit(n_q_points, number_of_values);
784 
785  handle = triangulation->register_data_attach(
786  std::bind(
788  this,
789  std::placeholders::_1,
790  std::placeholders::_2),
791  /*returns_variable_size_data=*/false);
792  }
793 
794 
795 
796  template <int dim, typename DataType>
797  void
799  {
800  triangulation->notify_ready_to_unpack(
801  handle,
802  std::bind(
804  this,
805  std::placeholders::_1,
806  std::placeholders::_2,
807  std::placeholders::_3));
808 
809  // invalidate the pointers
810  data_storage = nullptr;
811  triangulation = nullptr;
812  }
813 
814 
815 
816  template <int dim, typename DataType>
817  std::vector<char>
820  &cell,
822  dim>::CellStatus /*status*/)
823  {
824  pack_cell_data(cell, data_storage, matrix_quadrature);
825 
826  // project to FE
827  project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
828 
829  // to get consistent data sizes on each cell for the fixed size transfer,
830  // we won't allow compression
831  return Utilities::pack(matrix_dofs, /*allow_compression=*/false);
832  }
833 
834 
835 
836  template <int dim, typename DataType>
837  void
840  &cell,
841  const typename parallel::distributed::Triangulation<dim>::CellStatus
842  status,
843  const boost::iterator_range<std::vector<char>::const_iterator>
844  &data_range)
845  {
846  Assert((status !=
849  (void)status;
850 
851  matrix_dofs =
852  Utilities::unpack<FullMatrix<double>>(data_range.begin(),
853  data_range.end(),
854  /*allow_compression=*/false);
855 
856  if (cell->has_children())
857  {
858  // we need to first use prolongation matrix to get dofvalues on child
859  // cells based on dofvalues stored in the parent's data_store
860  for (unsigned int child = 0; child < cell->n_children(); ++child)
861  if (cell->child(child)->is_locally_owned())
862  {
863  projection_fe
864  ->get_prolongation_matrix(child, cell->refinement_case())
865  .mmult(matrix_dofs_child, matrix_dofs);
866 
867  // now we do the usual business of evaluating FE on quadrature
868  // points:
869  project_to_qp_matrix.mmult(matrix_quadrature,
870  matrix_dofs_child);
871 
872  // finally, put back into the map:
873  unpack_to_cell_data(cell->child(child),
874  matrix_quadrature,
875  data_storage);
876  }
877  }
878  else
879  {
880  // if there are no children, evaluate FE field at
881  // rhs_quadrature points.
882  project_to_qp_matrix.mmult(matrix_quadrature, matrix_dofs);
883 
884  // finally, put back into the map:
885  unpack_to_cell_data(cell, matrix_quadrature, data_storage);
886  }
887  }
888 
889  } // namespace distributed
890 
891 } // namespace parallel
892 
893 # endif // DEAL_II_WITH_P4EST
894 
895 #endif // DOXYGEN
896 DEAL_II_NAMESPACE_CLOSE
897 
898 #endif
size_type m() const
~CellDataStorage() override=default
void compute_interpolation_to_quadrature_points_matrix(const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, FullMatrix< double > &I_q)
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)
CellDataStorage()=default
STL namespace.
std::map< CellIteratorType, std::vector< std::shared_ptr< DataType > > > map
virtual unsigned int number_of_values() 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
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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)
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)
void prepare_for_coarsening_and_refinement(parallel::distributed::Triangulation< dim > &tria, CellDataStorage< CellIteratorType, DataType > &data_storage)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1174
virtual void unpack_values(const std::vector< double > &values)=0
virtual ~TransferableQuadraturePointData()=default
parallel::distributed::Triangulation< dim > * triangulation
std::vector< char > pack_function(const typename parallel::distributed::Triangulation< dim >::cell_iterator &cell, const typename parallel::distributed::Triangulation< dim >::CellStatus status)
static ::ExceptionBase & ExcNotImplemented()
typename parallel::distributed::Triangulation< dim >::cell_iterator CellIteratorType
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition: tria.h:269
virtual void pack_values(std::vector< double > &values) const =0
T max(const T &t, const MPI_Comm &mpi_communicator)
void initialize(const CellIteratorType &cell, const unsigned int number_of_data_points_per_cell)
static ::ExceptionBase & ExcInternalError()