16 #ifndef dealii__quadrature_point_data_h 17 #define dealii__quadrature_point_data_h 20 #include <deal.II/base/config.h> 22 #ifdef DEAL_II_WITH_CXX11 24 #include <deal.II/base/quadrature.h> 25 #include <deal.II/base/subscriptor.h> 26 #include <deal.II/grid/tria.h> 27 #include <deal.II/grid/tria_accessor.h> 28 #include <deal.II/distributed/tria.h> 29 #include <deal.II/fe/fe_tools.h> 30 #include <deal.II/fe/fe.h> 31 #include <deal.II/lac/vector.h> 35 #include <type_traits> 37 DEAL_II_NAMESPACE_OPEN
59 template <
typename CellIteratorType,
typename DataType>
94 template<
typename T=DataType>
96 const unsigned int number_of_data_points_per_cell);
103 template<
typename T=DataType>
104 void initialize(
const CellIteratorType &cell_start,
105 const CellIteratorType &cell_end,
106 const unsigned int number_of_data_points_per_cell);
117 bool erase(
const CellIteratorType &cell);
134 template<
typename T=DataType>
135 std::vector<std::shared_ptr<T> >
get_data(
const CellIteratorType &cell);
147 template<
typename T=DataType>
148 std::vector<std::shared_ptr<const T> >
get_data(
const CellIteratorType &cell)
const;
154 std::map<CellIteratorType, std::vector<std::shared_ptr<DataType>>>
map;
194 virtual unsigned int number_of_values()
const = 0;
205 virtual void pack_values(std::vector<double> &values)
const = 0;
214 virtual void unpack_values(
const std::vector<double> &values) = 0;
218 #ifdef DEAL_II_WITH_P4EST 221 namespace distributed
325 template <
int dim,
typename DataType>
329 static_assert(std::is_base_of<TransferableQuadraturePointData, DataType>::value,
330 "User's DataType class should be derived from QPData");
388 const typename parallel::distributed::Triangulation<dim,dim>::CellStatus status,
397 const typename parallel::distributed::Triangulation<dim,dim>::CellStatus status,
487 template <
typename CellIteratorType,
typename DataType>
496 template <
typename CellIteratorType,
typename DataType>
503 template <
typename CellIteratorType,
typename DataType>
506 const unsigned int n_q_points)
508 static_assert(std::is_base_of<DataType, T>::value,
509 "User's T class should be derived from user's DataType class");
511 if (
map.find(cell) ==
map.end())
513 map[cell] = std::vector<std::shared_ptr<DataType> >(n_q_points);
516 auto it =
map.find(cell);
517 for (
unsigned int q=0; q < n_q_points; q++)
518 it->second[q] = std::make_shared<T>();
524 template <
typename CellIteratorType,
typename DataType>
528 const unsigned int number)
531 if (it->is_locally_owned())
532 initialize<T>(it,number);
537 template <
typename CellIteratorType,
typename DataType>
540 const auto it =
map.find(cell);
541 for (
unsigned int i = 0; i < it->second.size(); i++)
543 Assert(it->second[i].unique(),
544 ExcMessage(
"Can not erase the cell data multiple objects reference its data."));
547 return (
map.erase(cell)==1);
553 template <
typename CellIteratorType,
typename DataType>
559 auto it =
map.begin();
560 while (it !=
map.end())
563 for (
unsigned int i = 0; i < it->second.size(); i++)
565 Assert(it->second[i].unique(),
566 ExcMessage(
"Can not erase the cell data, multiple objects reference it."));
575 template <
typename CellIteratorType,
typename DataType>
576 template <
typename T>
577 std::vector<std::shared_ptr<T> >
580 static_assert(std::is_base_of<DataType, T>::value,
581 "User's T class should be derived from user's DataType class");
583 auto it =
map.find(cell);
590 std::vector<std::shared_ptr<T>> res (it->second.size());
591 for (
unsigned int q = 0; q < res.size(); q++)
592 res[q] = std::dynamic_pointer_cast<T>(it->second[q]);
599 template <
typename CellIteratorType,
typename DataType>
600 template <
typename T>
601 std::vector<std::shared_ptr<const T> >
604 static_assert(std::is_base_of<DataType, T>::value,
605 "User's T class should be derived from user's DataType class");
607 auto it =
map.find(cell);
613 std::vector<std::shared_ptr<const T>> res (it->second.size());
614 for (
unsigned int q = 0; q < res.size(); q++)
615 res[q] = std::dynamic_pointer_cast<const T>(it->second[q]);
631 template<
typename CellIteratorType,
typename DataType>
637 static_assert(std::is_base_of<TransferableQuadraturePointData, DataType>::value,
638 "User's DataType class should be derived from QPData");
640 const std::vector<std::shared_ptr<const DataType> > qpd = data_storage->
get_data(cell);
642 const unsigned int n = matrix_data.
n();
644 std::vector<double> single_qp_data(n);
645 Assert (qpd.size() == matrix_data.
m(),
647 for (
unsigned int q = 0; q < qpd.size(); q++)
649 qpd[q]->pack_values(single_qp_data);
650 Assert (single_qp_data.size() == n,
653 for (
unsigned int i = 0; i < n; i++)
654 matrix_data(q,i) = single_qp_data[i];
664 template<
typename CellIteratorType,
typename DataType>
665 void unpack_to_cell_data
670 static_assert(std::is_base_of<TransferableQuadraturePointData, DataType>::value,
671 "User's DataType class should be derived from QPData");
673 std::vector<std::shared_ptr<DataType> > qpd = data_storage->
get_data(cell);
675 const unsigned int n = values_at_qp.
n();
677 std::vector<double> single_qp_data(n);
678 Assert(qpd.size() == values_at_qp.
m(),
681 for (
unsigned int q = 0; q < qpd.size(); q++)
683 for (
unsigned int i = 0; i < n; i++)
684 single_qp_data[i] = values_at_qp(q,i);
685 qpd[q]->unpack_values(single_qp_data);
690 #ifdef DEAL_II_WITH_P4EST 694 namespace distributed
696 template<
int dim,
typename DataType>
703 data_size_in_bytes(0),
704 n_q_points(rhs_quadrature.
size()),
705 project_to_fe_matrix(projection_fe->dofs_per_cell,n_q_points),
706 project_to_qp_matrix(n_q_points,projection_fe->dofs_per_cell),
708 data_storage(
nullptr),
709 triangulation(
nullptr)
711 Assert(projection_fe->n_components() == 1,
712 ExcMessage(
"ContinuousQuadratureDataTransfer requires scalar FiniteElement"));
715 (*projection_fe.get(),
718 project_to_fe_matrix);
721 (*projection_fe.get(),
723 project_to_qp_matrix);
728 template<
int dim,
typename DataType>
734 template<
int dim,
typename DataType>
739 Assert (data_storage ==
nullptr,
740 ExcMessage(
"This function can be called only once"));
741 triangulation = &tr_;
742 data_storage = &data_storage_;
744 unsigned int number_of_values = 0;
748 it != triangulation->end(); it++)
749 if (it->is_locally_owned())
751 std::vector<std::shared_ptr<DataType> > qpd = data_storage->
get_data(it);
752 number_of_values = qpd[0]->number_of_values();
756 number_of_values =
Utilities::MPI::max(number_of_values, triangulation->get_communicator ());
757 Assert (number_of_values > 0,
759 const unsigned int dofs_per_cell = projection_fe->dofs_per_cell;
760 matrix_dofs.reinit(dofs_per_cell,number_of_values);
761 matrix_dofs_child.reinit(dofs_per_cell,number_of_values);
762 matrix_quadrature.reinit(n_q_points,number_of_values);
763 data_size_in_bytes =
sizeof(double) * dofs_per_cell * number_of_values;
765 offset = triangulation->register_data_attach(data_size_in_bytes,
775 template<
int dim,
typename DataType>
778 triangulation->notify_ready_to_unpack(offset,
786 data_storage =
nullptr;
787 triangulation =
nullptr;
792 template<
int dim,
typename DataType>
795 const typename parallel::distributed::Triangulation<dim,dim>::CellStatus ,
798 double *data_store =
reinterpret_cast<double *
>(data);
800 pack_cell_data(cell,data_storage,matrix_quadrature);
803 project_to_fe_matrix.mmult(matrix_dofs, matrix_quadrature);
805 std::memcpy(data_store, &matrix_dofs(0,0), data_size_in_bytes);
810 template<
int dim,
typename DataType>
813 const typename parallel::distributed::Triangulation<dim,dim>::CellStatus status,
820 const double *data_store =
reinterpret_cast<const double *
>(data);
822 std::memcpy(&matrix_dofs(0,0), data_store, data_size_in_bytes);
824 if (cell->has_children())
828 for (
unsigned int child=0; child<cell->n_children(); ++child)
829 if (cell->child(child)->is_locally_owned())
831 projection_fe->get_prolongation_matrix(child, cell->refinement_case())
832 .mmult (matrix_dofs_child, matrix_dofs);
835 project_to_qp_matrix.mmult(matrix_quadrature,matrix_dofs_child);
838 unpack_to_cell_data(cell->child(child),matrix_quadrature,data_storage);
845 project_to_qp_matrix.mmult(matrix_quadrature,matrix_dofs);
848 unpack_to_cell_data(cell,matrix_quadrature,data_storage);
856 #endif // DEAL_II_WITH_P4EST 859 DEAL_II_NAMESPACE_CLOSE
std::map< CellIteratorType, std::vector< std::shared_ptr< DataType > > > map
std::vector< std::shared_ptr< T > > get_data(const CellIteratorType &cell)
bool erase(const CellIteratorType &cell)
unsigned int size() const
static::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const std::unique_ptr< const FiniteElement< dim > > projection_fe
std::size_t data_size_in_bytes
::Triangulation< dim, spacedim >::active_cell_iterator active_cell_iterator
FullMatrix< double > project_to_qp_matrix
parallel::distributed::Triangulation< dim >::cell_iterator CellIteratorType
virtual ~TransferableQuadraturePointData()
parallel::distributed::Triangulation< dim > * triangulation
TransferableQuadraturePointData()
::Triangulation< dim, spacedim >::cell_iterator cell_iterator
static::ExceptionBase & ExcNotImplemented()
CellDataStorage< CellIteratorType, DataType > * data_storage
FullMatrix< double > matrix_quadrature
FullMatrix< double > matrix_dofs_child
virtual FiniteElement< dim, spacedim > * clone() const =0
FullMatrix< double > project_to_fe_matrix
FullMatrix< double > matrix_dofs
T max(const T &t, const MPI_Comm &mpi_communicator)
void initialize(const CellIteratorType &cell, const unsigned int number_of_data_points_per_cell)
const unsigned int n_q_points
static::ExceptionBase & ExcInternalError()