44#include <boost/range/iterator_range_core.hpp>
63 template <
typename value_type>
66 const unsigned int dofs_per_cell)
68 for (
const auto &values : dof_values)
74 const std::size_t bytes_per_entry =
sizeof(value_type) * dofs_per_cell;
76 std::vector<char> buffer(dof_values.size() * bytes_per_entry);
77 for (
unsigned int i = 0; i < dof_values.size(); ++i)
78 std::memcpy(&buffer[i * bytes_per_entry],
90 template <
typename value_type>
91 std::vector<Vector<value_type>>
93 const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
94 const unsigned int dofs_per_cell)
96 const std::size_t bytes_per_entry =
sizeof(value_type) * dofs_per_cell;
97 const unsigned int n_elements = data_range.size() / bytes_per_entry;
101 std::vector<Vector<value_type>> unpacked_data;
102 unpacked_data.reserve(n_elements);
103 for (
unsigned int i = 0; i < n_elements; ++i)
106 std::memcpy(&dof_values(0),
107 &(*std::next(data_range.begin(), i * bytes_per_entry)),
109 unpacked_data.emplace_back(std::move(dof_values));
112 return unpacked_data;
118template <
int dim,
typename VectorType,
int spacedim>
121 const bool average_values)
122 : dof_handler(&dof, typeid(*this).name())
123 , average_values(average_values)
124 , handle(
numbers::invalid_unsigned_int)
129template <
int dim,
typename VectorType,
int spacedim>
133 const std::vector<const VectorType *> &all_in)
135 const ::internal::parallel::shared::
136 TemporarilyRestoreSubdomainIds<dim, spacedim>
137 subdomain_modifier(dof_handler->get_triangulation());
139 for (
unsigned int i = 0; i < all_in.size(); ++i)
140 Assert(all_in[i]->
size() == dof_handler->n_dofs(),
143 input_vectors = all_in;
144 register_data_attach();
149template <
int dim,
typename VectorType,
int spacedim>
154 std::vector<const VectorType *> temp(all_in.size());
156 for (std::size_t i = 0; i < temp.size(); ++i)
157 temp[i] = &(all_in[i]);
159 this->prepare_for_coarsening_and_refinement(temp);
164template <
int dim,
typename VectorType,
int spacedim>
170 &dof_handler->get_triangulation());
174 ExcMessage(
"You can only add one solution per "
175 "SolutionTransfer object."));
177 handle = tria->register_data_attach(
180 return this->pack_callback(cell_, status);
182 dof_handler->has_hp_capabilities());
187template <
int dim,
typename VectorType,
int spacedim>
192 std::vector<const VectorType *> all_in(1, &in);
193 prepare_for_coarsening_and_refinement(all_in);
198template <
int dim,
typename VectorType,
int spacedim>
201 const VectorType &in)
203 std::vector<const VectorType *> all_in(1, &in);
204 prepare_for_serialization(all_in);
209template <
int dim,
typename VectorType,
int spacedim>
212 const std::vector<const VectorType *> &all_in)
214 prepare_for_coarsening_and_refinement(all_in);
219template <
int dim,
typename VectorType,
int spacedim>
223 std::vector<VectorType *> all_in(1, &in);
229template <
int dim,
typename VectorType,
int spacedim>
232 std::vector<VectorType *> &all_in)
234 register_data_attach();
237 input_vectors.resize(all_in.size());
243template <
int dim,
typename VectorType,
int spacedim>
246 std::vector<VectorType *> &all_out)
248 const ::internal::parallel::shared::
249 TemporarilyRestoreSubdomainIds<dim, spacedim>
250 subdomain_modifier(dof_handler->get_triangulation());
252 Assert(input_vectors.size() == all_out.size(),
254 for (
unsigned int i = 0; i < all_out.size(); ++i)
255 Assert(all_out[i]->
size() == dof_handler->n_dofs(),
260 &dof_handler->get_triangulation());
265 "You can only call interpolate() once per SolutionTransfer object."));
267 using Number =
typename VectorType::value_type;
270 for (
auto *
const vec : all_out)
277 valence.reinit(*all_out[0]);
279 tria->notify_ready_to_unpack(
281 [
this, &all_out, &valence](
284 const boost::iterator_range<std::vector<char>::const_iterator>
286 this->unpack_callback(cell_, status, data_range, all_out, valence);
293 for (
const auto i : valence.locally_owned_elements())
294 valence[i] = (
static_cast<Number
>(valence[i]) == Number() ?
296 (Number(1.0) /
static_cast<Number
>(valence[i])));
299 for (
auto *
const vec : all_out)
308 for (
auto *
const vec : all_out)
312 input_vectors.clear();
317template <
int dim,
typename VectorType,
int spacedim>
320 std::vector<VectorType> &all_out)
322 std::vector<VectorType *> temp(all_out.size());
324 for (std::size_t i = 0; i < temp.size(); ++i)
325 temp[i] = &(all_out[i]);
327 this->interpolate(temp);
332template <
int dim,
typename VectorType,
int spacedim>
336 std::vector<VectorType *> all_out(1, &out);
337 interpolate(all_out);
342template <
int dim,
typename VectorType,
int spacedim>
351 std::vector<::Vector<typename VectorType::value_type>> dof_values(
352 input_vectors.size());
354 unsigned int fe_index = 0;
355 if (dof_handler->has_hp_capabilities())
362 fe_index = cell->future_fe_index();
372 for (
const auto &child : cell->child_iterators())
373 Assert(child->is_active() && child->coarsen_flag_set(),
374 typename ::Triangulation<
375 dim>::ExcInconsistentCoarseningFlags());
378 fe_index = ::internal::hp::DoFHandlerImplementation::
379 dominated_future_fe_on_children<dim, spacedim>(cell);
389 const unsigned int dofs_per_cell =
390 dof_handler->get_fe(fe_index).n_dofs_per_cell();
392 if (dofs_per_cell == 0)
393 return std::vector<char>();
395 auto it_input = input_vectors.cbegin();
396 auto it_output = dof_values.begin();
397 for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
399 it_output->reinit(dofs_per_cell);
400 cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
403 return pack_dof_values<typename VectorType::value_type>(dof_values,
409template <
int dim,
typename VectorType,
int spacedim>
414 const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
415 std::vector<VectorType *> &all_out,
420 unsigned int fe_index = 0;
421 if (dof_handler->has_hp_capabilities())
428 fe_index = cell->active_fe_index();
439 fe_index = cell->child(0)->active_fe_index();
440 for (
unsigned int child_index = 1;
441 child_index < cell->n_children();
443 Assert(cell->child(child_index)->active_fe_index() == fe_index,
454 const unsigned int dofs_per_cell =
455 dof_handler->get_fe(fe_index).n_dofs_per_cell();
457 if (dofs_per_cell == 0)
460 const std::vector<::Vector<typename VectorType::value_type>>
462 unpack_dof_values<typename VectorType::value_type>(data_range,
470 for (
auto it_dof_values = dof_values.begin();
471 it_dof_values != dof_values.end();
474 dofs_per_cell == it_dof_values->size(),
476 "The transferred data was packed with a different number of dofs than the "
477 "currently registered FE object assigned to the DoFHandler has."));
480 auto it_input = dof_values.cbegin();
481 auto it_output = all_out.begin();
482 for (; it_input != dof_values.cend(); ++it_input, ++it_output)
484 cell->distribute_local_to_global_by_interpolation(*it_input,
488 cell->set_dof_values_by_interpolation(*it_input,
498 cell->distribute_local_to_global_by_interpolation(ones,
506template <
int dim,
typename VectorType,
int spacedim>
518 template <
int dim,
typename VectorType,
int spacedim>
521 : dof_handler(&dof, typeid(*this).name())
529 "You are calling the ::SolutionTransfer class "
530 "with a DoFHandler that is built on a "
531 "parallel::distributed::Triangulation. This will not "
532 "work for parallel computations. You probably want to "
533 "use the parallel::distributed::SolutionTransfer class."));
538 template <
int dim,
typename VectorType,
int spacedim>
546 template <
int dim,
typename VectorType,
int spacedim>
550 indices_on_cell.clear();
551 dof_values_on_cell.clear();
559 template <
int dim,
typename VectorType,
int spacedim>
563 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
564 Assert(prepared_for != coarsening_and_refinement,
565 ExcAlreadyPrepForCoarseAndRef());
577 const ::internal::parallel::shared::
578 TemporarilyRestoreSubdomainIds<dim, spacedim>
579 subdomain_modifier(dof_handler->get_triangulation());
581 const unsigned int n_active_cells =
582 dof_handler->get_triangulation().n_active_cells();
583 n_dofs_old = dof_handler->n_dofs();
586 std::vector<std::vector<types::global_dof_index>>(n_active_cells)
587 .
swap(indices_on_cell);
589 for (
const auto &cell : dof_handler->active_cell_iterators())
591 const unsigned int i = cell->active_cell_index();
592 indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
599 cell->get_dof_indices(indices_on_cell[i]);
600 cell_map[std::make_pair(cell->level(), cell->index())] =
603 prepared_for = pure_refinement;
608 template <
int dim,
typename VectorType,
int spacedim>
611 const VectorType &in,
612 VectorType &out)
const
614 Assert(prepared_for == pure_refinement, ExcNotPrepared());
615 Assert(in.size() == n_dofs_old,
617 Assert(out.size() == dof_handler->n_dofs(),
620 ExcMessage(
"Vectors cannot be used as input and output"
621 " at the same time!"));
631 const ::internal::parallel::shared::
632 TemporarilyRestoreSubdomainIds<dim, spacedim>
633 subdomain_modifier(dof_handler->get_triangulation());
637 typename std::map<std::pair<unsigned int, unsigned int>,
639 cell_map_end = cell_map.end();
641 for (
const auto &cell : dof_handler->cell_iterators())
644 cell_map.find(std::make_pair(cell->level(), cell->index()));
646 if (pointerstruct != cell_map_end)
654 const unsigned int this_fe_index =
656 const unsigned int dofs_per_cell =
657 cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
658 local_values.
reinit(dofs_per_cell,
true);
663 Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
665 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
668 in, (*pointerstruct->second.indices_ptr)[i]);
669 cell->set_dof_values_by_interpolation(local_values,
694 template <
int dim,
int spacedim>
703 const ::hp::FECollection<dim, spacedim> &fe =
705 matrices.reinit(fe.size(), fe.size());
706 for (
unsigned int i = 0; i < fe.size(); ++i)
707 for (
unsigned int j = 0; j < fe.size(); ++j)
710 matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
711 fe[j].n_dofs_per_cell());
722 fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
725 ExcInterpolationNotImplemented &)
727 matrices(i, j).reinit(0, 0);
733 template <
int dim,
int spacedim>
736 std::vector<std::vector<bool>> &)
739 template <
int dim,
int spacedim>
742 const ::hp::FECollection<dim, spacedim> &fe,
743 std::vector<std::vector<bool>> &restriction_is_additive)
745 restriction_is_additive.resize(fe.size());
746 for (
unsigned int f = 0; f < fe.size(); ++f)
748 restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
749 for (
unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
750 restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
757 template <
int dim,
typename VectorType,
int spacedim>
762 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
763 Assert(prepared_for != coarsening_and_refinement,
764 ExcAlreadyPrepForCoarseAndRef());
767 n_dofs_old = dof_handler->n_dofs();
768 const unsigned int in_size = all_in.size();
772 ExcMessage(
"The array of input vectors you pass to this "
773 "function has no elements. This is not useful."));
774 for (
unsigned int i = 0; i < in_size; ++i)
789 const ::internal::parallel::shared::
790 TemporarilyRestoreSubdomainIds<dim, spacedim>
791 subdomain_modifier(dof_handler->get_triangulation());
796 unsigned int n_cells_to_coarsen = 0;
797 unsigned int n_cells_to_stay_or_refine = 0;
798 for (
const auto &act_cell : dof_handler->active_cell_iterators())
800 if (act_cell->coarsen_flag_set())
801 ++n_cells_to_coarsen;
803 ++n_cells_to_stay_or_refine;
805 Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
806 dof_handler->get_triangulation().n_active_cells(),
809 unsigned int n_coarsen_fathers = 0;
810 for (
const auto &cell : dof_handler->cell_iterators())
811 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
814 (void)n_cells_to_coarsen;
819 std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
820 .
swap(indices_on_cell);
822 std::vector<std::vector<Vector<typename VectorType::value_type>>>(
824 std::vector<Vector<typename VectorType::value_type>>(in_size))
825 .swap(dof_values_on_cell);
828 std::vector<std::vector<bool>> restriction_is_additive;
832 restriction_is_additive);
837 unsigned int n_sr = 0, n_cf = 0;
838 for (
const auto &cell : dof_handler->cell_iterators())
841 if (cell->is_active() && !cell->coarsen_flag_set())
843 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
844 indices_on_cell[n_sr].resize(dofs_per_cell);
849 cell->get_dof_indices(indices_on_cell[n_sr]);
850 cell_map[std::make_pair(cell->level(), cell->index())] =
851 Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
856 else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
864 std::set<unsigned int> fe_indices_children;
865 for (
const auto &child : cell->child_iterators())
867 Assert(child->is_active() && child->coarsen_flag_set(),
868 typename ::Triangulation<
869 dim>::ExcInconsistentCoarseningFlags());
871 fe_indices_children.insert(child->active_fe_index());
875 const unsigned int target_fe_index =
876 dof_handler->get_fe_collection().find_dominated_fe_extended(
877 fe_indices_children, 0);
880 ::internal::hp::DoFHandlerImplementation::
881 ExcNoDominatedFiniteElementOnChildren());
883 const unsigned int dofs_per_cell =
884 dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
886 std::vector<Vector<typename VectorType::value_type>>(
888 .
swap(dof_values_on_cell[n_cf]);
896 for (
unsigned int j = 0; j < in_size; ++j)
897 cell->get_interpolated_dof_values(all_in[j],
898 dof_values_on_cell[n_cf][j],
900 cell_map[std::make_pair(cell->level(), cell->index())] =
908 prepared_for = coarsening_and_refinement;
913 template <
int dim,
typename VectorType,
int spacedim>
918 std::vector<VectorType> all_in(1, in);
919 prepare_for_coarsening_and_refinement(all_in);
924 template <
int dim,
typename VectorType,
int spacedim>
927 const std::vector<VectorType> &all_in,
928 std::vector<VectorType> &all_out)
const
930 const unsigned int size = all_in.size();
932 Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
934 for (
unsigned int i = 0; i <
size; ++i)
937 for (
unsigned int i = 0; i < all_out.size(); ++i)
938 Assert(all_out[i].
size() == dof_handler->n_dofs(),
940 for (
unsigned int i = 0; i <
size; ++i)
941 for (
unsigned int j = 0; j <
size; ++j)
942 Assert(&all_in[i] != &all_out[j],
943 ExcMessage(
"Vectors cannot be used as input and output"
944 " at the same time!"));
955 const ::internal::parallel::shared::
956 TemporarilyRestoreSubdomainIds<dim, spacedim>
957 subdomain_modifier(dof_handler->get_triangulation());
960 std::vector<types::global_dof_index> dofs;
962 typename std::map<std::pair<unsigned int, unsigned int>,
964 cell_map_end = cell_map.end();
970 for (
const auto &cell : dof_handler->cell_iterators())
973 cell_map.find(std::make_pair(cell->level(), cell->index()));
975 if (pointerstruct != cell_map_end)
977 const std::vector<types::global_dof_index> *
const indexptr =
978 pointerstruct->second.indices_ptr;
980 const std::vector<Vector<typename VectorType::value_type>>
981 *
const valuesptr = pointerstruct->second.dof_values_ptr;
984 if (indexptr !=
nullptr)
988 const unsigned int old_fe_index =
989 pointerstruct->second.active_fe_index;
993 unsigned int in_size = indexptr->size();
994 for (
unsigned int j = 0; j <
size; ++j)
996 tmp.
reinit(in_size,
true);
997 for (
unsigned int i = 0; i < in_size; ++i)
999 all_in[j], (*indexptr)[i]);
1001 cell->set_dof_values_by_interpolation(tmp,
1013 const unsigned int dofs_per_cell =
1014 cell->get_fe().n_dofs_per_cell();
1015 dofs.resize(dofs_per_cell);
1018 cell->get_dof_indices(dofs);
1021 for (
unsigned int j = 0; j <
size; ++j)
1029 const unsigned int active_fe_index =
1030 cell->active_fe_index();
1031 if (active_fe_index !=
1032 pointerstruct->second.active_fe_index)
1034 const unsigned int old_index =
1035 pointerstruct->second.active_fe_index;
1037 interpolation_hp(active_fe_index, old_index);
1040 if (interpolation_matrix.empty())
1041 tmp.
reinit(dofs_per_cell,
false);
1044 tmp.
reinit(dofs_per_cell,
true);
1046 interpolation_matrix.
n());
1048 interpolation_matrix.
m());
1049 interpolation_matrix.
vmult(tmp, (*valuesptr)[j]);
1054 data = &(*valuesptr)[j];
1057 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1059 (*
data)(i), dofs[i], all_out[j]);
1070 for (
auto &vec : all_out)
1076 template <
int dim,
typename VectorType,
int spacedim>
1079 const VectorType &in,
1080 VectorType &out)
const
1082 Assert(in.size() == n_dofs_old,
1084 Assert(out.size() == dof_handler->n_dofs(),
1087 std::vector<VectorType> all_in = {in};
1088 std::vector<VectorType> all_out = {out};
1090 interpolate(all_in, all_out);
1097 template <
int dim,
typename VectorType,
int spacedim>
1107 sizeof(prepared_for) +
1114 template <
int dim,
typename VectorType,
int spacedim>
1119 return sizeof(*this);
1126#define SPLIT_INSTANTIATIONS_COUNT 4
1127#ifndef SPLIT_INSTANTIATIONS_INDEX
1128# define SPLIT_INSTANTIATIONS_INDEX 0
1130#include "numerics/solution_transfer.inst"
@ children_will_be_coarsened
const hp::FECollection< dim, spacedim > & get_fe_collection() const
bool has_hp_capabilities() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void refine_interpolate(const VectorType &in, VectorType &out) const
void interpolate(const std::vector< VectorType > &all_in, std::vector< VectorType > &all_out) const
std::size_t memory_consumption() const
ObserverPointer< const DoFHandler< dim, spacedim >, SolutionTransfer< dim, VectorType, spacedim > > dof_handler
void prepare_for_coarsening_and_refinement(const std::vector< VectorType > &all_in)
void prepare_for_pure_refinement()
SolutionTransfer(const DoFHandler< dim, spacedim > &dof)
std::vector< char > pack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellStatus status)
void unpack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellStatus status, const boost::iterator_range< std::vector< char >::const_iterator > &data_range, std::vector< VectorType * > &all_out, VectorType &valence)
void register_data_attach()
void deserialize(VectorType &in)
void prepare_for_coarsening_and_refinement(const std::vector< const VectorType * > &all_in)
void prepare_for_serialization(const VectorType &in)
void interpolate(std::vector< VectorType * > &all_out)
SolutionTransfer(const DoFHandler< dim, spacedim > &dof_handler, const bool average_values=false)
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
std::vector< index_type > data
void restriction_additive(const FiniteElement< dim, spacedim > &, std::vector< std::vector< bool > > &)
void extract_interpolation_matrices(const DoFHandler< dim, spacedim > &dof, ::Table< 2, FullMatrix< double > > &matrices)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
constexpr unsigned int invalid_unsigned_int
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
std::size_t memory_consumption() const
unsigned int active_fe_index
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)