17 #include <deal.II/distributed/tria.h> 19 #include <deal.II/dofs/dof_tools.h> 21 #include <deal.II/fe/fe_tools.h> 23 #include <deal.II/matrix_free/shape_info.h> 25 #include <deal.II/multigrid/mg_transfer_internal.h> 27 DEAL_II_NAMESPACE_OPEN
42 DoFPair(
const unsigned int level,
47 , level_dof_index(level_dof_index)
51 : level(
numbers::invalid_unsigned_int)
59 template <
int dim,
int spacedim>
62 const ::DoFHandler<dim, spacedim> &mg_dof,
64 std::vector<std::vector<
65 std::pair<types::global_dof_index, types::global_dof_index>>>
67 std::vector<std::vector<
68 std::pair<types::global_dof_index, types::global_dof_index>>>
69 ©_indices_global_mine,
70 std::vector<std::vector<
71 std::pair<types::global_dof_index, types::global_dof_index>>>
72 & copy_indices_level_mine,
73 const bool skip_interface_dofs)
86 std::vector<DoFPair> send_data_temp;
88 const unsigned int n_levels =
89 mg_dof.get_triangulation().n_global_levels();
90 copy_indices.resize(n_levels);
91 copy_indices_global_mine.resize(n_levels);
92 copy_indices_level_mine.resize(n_levels);
96 const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
97 std::vector<types::global_dof_index> global_dof_indices(dofs_per_cell);
98 std::vector<types::global_dof_index> level_dof_indices(dofs_per_cell);
100 for (
unsigned int level = 0; level < n_levels; ++level)
102 std::vector<bool> dof_touched(globally_relevant.
n_elements(),
false);
103 copy_indices[level].clear();
104 copy_indices_level_mine[level].clear();
105 copy_indices_global_mine[level].clear();
107 typename ::DoFHandler<dim, spacedim>::active_cell_iterator
108 level_cell = mg_dof.begin_active(level);
109 const typename ::DoFHandler<dim, spacedim>::active_cell_iterator
110 level_end = mg_dof.end_active(level);
112 for (; level_cell != level_end; ++level_cell)
114 if (mg_dof.get_triangulation().locally_owned_subdomain() !=
116 (level_cell->level_subdomain_id() ==
118 level_cell->subdomain_id() ==
124 level_cell->get_dof_indices(global_dof_indices);
125 level_cell->get_mg_dof_indices(level_dof_indices);
127 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
131 if (skip_interface_dofs && mg_constrained_dofs !=
nullptr &&
133 level, level_dof_indices[i]))
140 if (dof_touched[global_idx])
142 bool global_mine = mg_dof.locally_owned_dofs().is_element(
143 global_dof_indices[i]);
145 mg_dof.locally_owned_mg_dofs(level).is_element(
146 level_dof_indices[i]);
149 if (global_mine && level_mine)
151 copy_indices[level].emplace_back(global_dof_indices[i],
152 level_dof_indices[i]);
154 else if (global_mine)
156 copy_indices_global_mine[level].emplace_back(
157 global_dof_indices[i], level_dof_indices[i]);
160 send_data_temp.emplace_back(level,
161 global_dof_indices[i],
162 level_dof_indices[i]);
169 dof_touched[global_idx] =
true;
174 const ::parallel::Triangulation<dim, spacedim> *tria =
176 &mg_dof.get_triangulation()));
178 send_data_temp.size() == 0 || tria !=
nullptr,
180 "We should only be sending information with a parallel Triangulation!"));
182 #ifdef DEAL_II_WITH_MPI 190 const std::set<types::subdomain_id> &neighbors =
191 tria->level_ghost_owners();
192 std::map<int, std::vector<DoFPair>> send_data;
196 for (
const auto &dofpair : send_data_temp)
198 std::set<types::subdomain_id>::iterator it;
199 for (it = neighbors.begin(); it != neighbors.end(); ++it)
202 .locally_owned_mg_dofs_per_processor(dofpair.level)[*it]
203 .is_element(dofpair.level_dof_index))
205 send_data[*it].push_back(dofpair);
211 Assert(it != neighbors.end(),
216 std::vector<MPI_Request> requests;
218 for (
const auto dest : neighbors)
220 requests.push_back(MPI_Request());
221 std::vector<DoFPair> &data = send_data[dest];
227 const int ierr = MPI_Isend(data.data(),
228 data.size() *
sizeof(data[0]),
232 tria->get_communicator(),
233 &*requests.rbegin());
238 const int ierr = MPI_Isend(
nullptr,
243 tria->get_communicator(),
244 &*requests.rbegin());
253 std::vector<DoFPair> receive_buffer;
254 for (
unsigned int counter = 0; counter < neighbors.size();
259 int ierr = MPI_Probe(MPI_ANY_SOURCE,
261 tria->get_communicator(),
264 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
269 ierr = MPI_Recv(
nullptr,
274 tria->get_communicator(),
280 int count = len /
sizeof(DoFPair);
281 Assert(static_cast<int>(count *
sizeof(DoFPair)) == len,
283 receive_buffer.resize(count);
285 void *ptr = receive_buffer.data();
291 tria->get_communicator(),
295 for (
const auto &dof_pair : receive_buffer)
297 copy_indices_level_mine[dof_pair.level].emplace_back(
298 dof_pair.global_dof_index, dof_pair.level_dof_index);
304 if (requests.size() > 0)
306 const int ierr = MPI_Waitall(requests.size(),
308 MPI_STATUSES_IGNORE);
316 const int ierr = MPI_Barrier(tria->get_communicator());
325 std::less<std::pair<types::global_dof_index, types::global_dof_index>>
327 for (
auto &level_indices : copy_indices)
328 std::sort(level_indices.begin(), level_indices.end(), compare);
329 for (
auto &level_indices : copy_indices_level_mine)
330 std::sort(level_indices.begin(), level_indices.end(), compare);
331 for (
auto &level_indices : copy_indices_global_mine)
332 std::sort(level_indices.begin(), level_indices.end(), compare);
339 template <
typename Number>
341 reinit_ghosted_vector(
343 std::vector<types::global_dof_index> & ghosted_level_dofs,
344 const MPI_Comm & communicator,
346 std::vector<std::pair<unsigned int, unsigned int>>
347 ©_indices_global_mine)
349 std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
351 ghosted_dofs.
add_indices(ghosted_level_dofs.begin(),
352 std::unique(ghosted_level_dofs.begin(),
353 ghosted_level_dofs.end()));
354 ghosted_dofs.compress();
357 if (ghosted_level_vector.
size() == locally_owned.
size())
362 ghosted_dofs.add_indices(part->ghost_indices());
363 for (
auto &indices : copy_indices_global_mine)
365 ghosted_dofs.index_within_set(
366 part->local_to_global(indices.second));
368 ghosted_level_vector.
reinit(locally_owned, ghosted_dofs, communicator);
373 copy_indices_to_mpi_local_numbers(
375 const std::vector<types::global_dof_index> &mine,
376 const std::vector<types::global_dof_index> &remote,
377 std::vector<unsigned int> & localized_indices)
379 localized_indices.resize(mine.size() + remote.size(),
381 for (
unsigned int i = 0; i < mine.size(); ++i)
385 for (
unsigned int i = 0; i < remote.size(); ++i)
394 compute_shift_within_children(
const unsigned int child,
395 const unsigned int fe_shift_1d,
396 const unsigned int fe_degree)
400 unsigned int c_tensor_index[dim];
401 unsigned int tmp = child;
402 for (
unsigned int d = 0;
d < dim; ++
d)
404 c_tensor_index[
d] = tmp % 2;
407 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
408 unsigned int factor = 1;
409 unsigned int shift = fe_shift_1d * c_tensor_index[0];
410 for (
unsigned int d = 1;
d < dim; ++
d)
412 factor *= n_child_dofs_1d;
413 shift = shift + factor * fe_shift_1d * c_tensor_index[
d];
423 const unsigned int child,
424 const unsigned int fe_shift_1d,
425 const unsigned int fe_degree,
426 const std::vector<unsigned int> & lexicographic_numbering,
427 const std::vector<types::global_dof_index> &local_dof_indices,
430 const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
431 const unsigned int shift =
432 compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
434 local_dof_indices.size() / Utilities::fixed_power<dim>(fe_degree + 1);
436 const unsigned int n_scalar_cell_dofs =
437 Utilities::fixed_power<dim>(n_child_dofs_1d);
439 for (
unsigned int k = 0; k < (dim > 2 ? (fe_degree + 1) : 1); ++k)
440 for (
unsigned int j = 0; j < (dim > 1 ? (fe_degree + 1) : 1); ++j)
441 for (
unsigned int i = 0; i < (fe_degree + 1); ++i, ++m)
443 const unsigned int index =
444 c * n_scalar_cell_dofs +
445 k * n_child_dofs_1d * n_child_dofs_1d + j * n_child_dofs_1d +
449 local_dof_indices[lexicographic_numbering[m]],
451 indices[index] = local_dof_indices[lexicographic_numbering[m]];
455 template <
int dim,
typename Number>
457 setup_element_info(ElementInfo<Number> & elem_info,
459 const ::DoFHandler<dim> &mg_dof)
462 elem_info.n_components = mg_dof.get_fe().element_multiplicity(0);
464 elem_info.n_components,
465 mg_dof.get_fe().dofs_per_cell);
467 elem_info.fe_degree = fe.
degree;
492 const unsigned int n_child_dofs_1d =
496 elem_info.n_child_cell_dofs =
497 elem_info.n_components * Utilities::fixed_power<dim>(n_child_dofs_1d);
501 shape_info.
reinit(dummy_quadrature, mg_dof.get_fe(), 0);
508 for (
unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
512 .prolongation_matrix_1d[i * n_child_dofs_1d + j + c * shift] =
525 const unsigned int level,
526 std::vector<types::global_dof_index> &dof_indices)
528 if (mg_constrained_dofs !=
nullptr &&
531 for (
auto &ind : dof_indices)
549 template <
int dim,
typename Number>
552 const ::DoFHandler<dim> & mg_dof,
554 ElementInfo<Number> & elem_info,
555 std::vector<std::vector<unsigned int>> &level_dof_indices,
556 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
557 & parent_child_connect,
558 std::vector<unsigned int> &n_owned_level_cells,
559 std::vector<std::vector<std::vector<unsigned short>>> &dirichlet_indices,
560 std::vector<std::vector<Number>> & weights_on_refined,
561 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
562 ©_indices_global_mine,
564 &ghosted_level_vector)
566 level_dof_indices.clear();
567 parent_child_connect.clear();
568 n_owned_level_cells.clear();
569 dirichlet_indices.clear();
570 weights_on_refined.clear();
576 const ::Triangulation<dim> &tria = mg_dof.get_triangulation();
582 std::string fe_name = mg_dof.get_fe().base_element(0).get_name();
584 const std::size_t template_starts = fe_name.find_first_of(
'<');
585 Assert(fe_name[template_starts + 1] ==
586 (dim == 1 ?
'1' : (dim == 2 ?
'2' :
'3')),
588 fe_name[template_starts + 1] =
'1';
590 const std::unique_ptr<FiniteElement<1>> fe(
591 FETools::get_fe_by_name<1, 1>(fe_name));
593 setup_element_info(elem_info, *fe, mg_dof);
598 const unsigned int n_levels = tria.n_global_levels();
599 level_dof_indices.resize(n_levels);
600 parent_child_connect.resize(n_levels - 1);
601 n_owned_level_cells.resize(n_levels - 1);
602 std::vector<std::vector<unsigned int>> coarse_level_indices(n_levels - 1);
603 for (
unsigned int level = 0;
604 level < std::min(tria.n_levels(), n_levels - 1);
606 coarse_level_indices[level].resize(tria.n_raw_cells(level),
608 std::vector<types::global_dof_index> local_dof_indices(
609 mg_dof.get_fe().dofs_per_cell);
610 dirichlet_indices.resize(n_levels - 1);
615 if (ghosted_level_vector.max_level() != n_levels - 1)
616 ghosted_level_vector.resize(0, n_levels - 1);
618 for (
unsigned int level = n_levels - 1; level > 0; --level)
620 unsigned int counter = 0;
621 std::vector<types::global_dof_index> global_level_dof_indices;
622 std::vector<types::global_dof_index> global_level_dof_indices_remote;
623 std::vector<types::global_dof_index> ghosted_level_dofs;
624 std::vector<types::global_dof_index> global_level_dof_indices_l0;
625 std::vector<types::global_dof_index> ghosted_level_dofs_l0;
628 typename ::DoFHandler<dim>::cell_iterator cell,
629 endc = mg_dof.end(level - 1);
630 for (cell = mg_dof.begin(level - 1); cell != endc; ++cell)
634 if (!cell->has_children())
637 bool consider_cell =
false;
638 if (tria.locally_owned_subdomain() ==
640 cell->level_subdomain_id() == tria.locally_owned_subdomain())
641 consider_cell =
true;
646 bool cell_is_remote = !consider_cell;
647 for (
unsigned int c = 0;
648 c < GeometryInfo<dim>::max_children_per_cell;
650 if (cell->child(c)->level_subdomain_id() ==
651 tria.locally_owned_subdomain())
653 consider_cell =
true;
668 std::vector<types::global_dof_index> &next_indices =
669 cell_is_remote ? global_level_dof_indices_remote :
670 global_level_dof_indices;
671 const std::size_t start_index = next_indices.size();
672 next_indices.resize(start_index + elem_info.n_child_cell_dofs,
674 for (
unsigned int c = 0;
675 c < GeometryInfo<dim>::max_children_per_cell;
678 if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
679 tria.locally_owned_subdomain())
681 cell->child(c)->get_mg_dof_indices(local_dof_indices);
683 replace(mg_constrained_dofs, level, local_dof_indices);
686 mg_dof.locally_owned_mg_dofs(level);
687 for (
const auto local_dof_index : local_dof_indices)
688 if (!owned_level_dofs.
is_element(local_dof_index))
689 ghosted_level_dofs.push_back(local_dof_index);
691 add_child_indices<dim>(c,
695 elem_info.lexicographic_numbering,
697 &next_indices[start_index]);
700 if (cell->child(c)->has_children() &&
701 (tria.locally_owned_subdomain() ==
703 cell->child(c)->level_subdomain_id() ==
704 tria.locally_owned_subdomain()))
706 const unsigned int child_index =
707 coarse_level_indices[level][cell->child(c)->index()];
709 parent_child_connect[level].size());
710 unsigned int parent_index = counter;
719 start_index / elem_info.n_child_cell_dofs +
721 parent_child_connect[level][child_index] =
722 std::make_pair(parent_index, c);
724 static_cast<unsigned short>(-1));
728 if (mg_constrained_dofs !=
nullptr)
729 for (
unsigned int i = 0;
730 i < mg_dof.get_fe().dofs_per_cell;
735 [elem_info.lexicographic_numbering[i]]))
736 dirichlet_indices[level][child_index].push_back(i);
742 coarse_level_indices[level - 1].size());
743 coarse_level_indices[level - 1][cell->index()] = counter++;
750 if (level == 1 && !cell_is_remote)
752 cell->get_mg_dof_indices(local_dof_indices);
754 replace(mg_constrained_dofs, level - 1, local_dof_indices);
756 const IndexSet &owned_level_dofs_l0 =
757 mg_dof.locally_owned_mg_dofs(0);
758 for (
const auto local_dof_index : local_dof_indices)
759 if (!owned_level_dofs_l0.
is_element(local_dof_index))
760 ghosted_level_dofs_l0.push_back(local_dof_index);
762 const std::size_t start_index =
763 global_level_dof_indices_l0.size();
764 global_level_dof_indices_l0.resize(
765 start_index + elem_info.n_child_cell_dofs,
767 add_child_indices<dim>(
771 elem_info.lexicographic_numbering,
773 &global_level_dof_indices_l0[start_index]);
775 dirichlet_indices[0].emplace_back();
776 if (mg_constrained_dofs !=
nullptr)
777 for (
unsigned int i = 0; i < mg_dof.get_fe().dofs_per_cell;
781 local_dof_indices[elem_info
782 .lexicographic_numbering[i]]))
783 dirichlet_indices[0].back().push_back(i);
791 global_level_dof_indices.size());
792 n_owned_level_cells[level - 1] = counter;
793 dirichlet_indices[level - 1].resize(counter);
794 parent_child_connect[level - 1].resize(
802 if (level < n_levels - 1)
803 for (std::vector<std::pair<unsigned int, unsigned int>>::iterator
804 i = parent_child_connect[level].begin();
805 i != parent_child_connect[level].end();
807 if (i->first >= tria.n_cells(level))
809 i->first -= tria.n_cells(level);
816 const MPI_Comm communicator =
819 reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(level),
822 ghosted_level_vector[level],
823 copy_indices_global_mine[level]);
825 copy_indices_to_mpi_local_numbers(
826 *ghosted_level_vector[level].get_partitioner(),
827 global_level_dof_indices,
828 global_level_dof_indices_remote,
829 level_dof_indices[level]);
833 for (
unsigned int i = 0; i < parent_child_connect[0].size(); ++i)
834 parent_child_connect[0][i] = std::make_pair(i, 0U);
836 reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(0),
837 ghosted_level_dofs_l0,
839 ghosted_level_vector[0],
840 copy_indices_global_mine[0]);
842 copy_indices_to_mpi_local_numbers(
843 *ghosted_level_vector[0].get_partitioner(),
844 global_level_dof_indices_l0,
845 std::vector<types::global_dof_index>(),
846 level_dof_indices[0]);
852 const unsigned int n_child_dofs_1d =
857 weights_on_refined.resize(n_levels - 1);
858 for (
unsigned int level = 1; level < n_levels; ++level)
860 ghosted_level_vector[level] = 0;
861 for (
unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
862 for (
unsigned int j = 0; j < elem_info.n_child_cell_dofs; ++j)
863 ghosted_level_vector[level].local_element(
864 level_dof_indices[level][elem_info.n_child_cell_dofs * c +
867 ghosted_level_vector[level].update_ghost_values();
869 std::vector<unsigned int> degree_to_3(n_child_dofs_1d);
871 for (
unsigned int i = 1; i < n_child_dofs_1d - 1; ++i)
873 degree_to_3.back() = 2;
877 weights_on_refined[level - 1].resize(n_owned_level_cells[level - 1] *
878 Utilities::fixed_power<dim>(3));
879 for (
unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
880 for (
unsigned int k = 0, m = 0; k < (dim > 2 ? n_child_dofs_1d : 1);
882 for (
unsigned int j = 0; j < (dim > 1 ? n_child_dofs_1d : 1); ++j)
884 unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j];
885 for (
unsigned int i = 0; i < n_child_dofs_1d; ++i, ++m)
886 weights_on_refined[level -
887 1][c * Utilities::fixed_power<dim>(3) +
888 shift + degree_to_3[i]] =
890 ghosted_level_vector[level].local_element(
891 level_dof_indices[level]
892 [elem_info.n_child_cell_dofs * c + m]);
902 #include "mg_transfer_internal.inst" 904 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
#define AssertDimension(dim1, dim2)
const types::subdomain_id invalid_subdomain_id
void reinit(const Quadrature< 1 > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
#define AssertIndexRange(index, range)
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
const unsigned int degree
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
#define AssertThrow(cond, exc)
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
unsigned int global_to_local(const types::global_dof_index global_index) const
const unsigned int dofs_per_line
size_type n_constraints() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
size_type index_within_set(const size_type global_index) const
virtual MPI_Comm get_communicator() const
bool is_identity_constrained(const size_type line_n) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const unsigned int dofs_per_cell
const std::vector< Point< dim > > & get_unit_support_points() const
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
unsigned int global_dof_index
const types::subdomain_id artificial_subdomain_id
virtual size_type size() const override
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
#define AssertThrowMPI(error_code)
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
static ::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
bool is_element(const size_type index) const
std::vector< unsigned int > lexicographic_numbering
const types::global_dof_index invalid_dof_index
size_type n_elements() const
static ::ExceptionBase & ExcInternalError()